
# make sure this is the directory where the files are:
dir1 <- "C:\\Michael\\Research\\rep\\"

# replication options (set to false unless you want to wait a long time)
fullrep_ideal <- FALSE

# load an r library, downloading if necesary
library0 <- function(s)
{
	if(s %in% rownames(installed.packages()) == FALSE) {install.packages(s,repos='http://cran.us.r-project.org')}
	require(s,character.only=TRUE)
}

# install michael ideal point estimation package
#install.packages("https://sites.google.com/a/stonybrook.edu/mperess/ipe_PSRM.tar.gz")

library0("ipe")
library0("optrees")
library0("openxlsx") # for writing excel files
library0("fastmatch")
library0("stringr")
library0("Formula") # for formulas with multiple parts
library0("readxl")
library0("numDeriv") # for numerical derivatives
library0("MASS") # for ordered probit

# some useful functions
check.integer <- function(x) { x == round(x) }

print0 <- function(s1,s2=logical(0))
{
	if(length(s2) == 0)
	{
		print(s1)
		flush.console()
	} else {
		print(paste(s1,": ",s2,sep=""))
		flush.console()
	}
}

apply0 <- function(func1,x,subset=numeric(0),by=numeric(0),bys=numeric(0),...)
{
	if(length(subset) == 0) {
		x1 <- subset(x,!is.na(x))
		if(length(by) != 0) by1 <- subset(by,!is.na(x))
	} else if(length(subset) != length(x)) {
		# Error
		stop("Incompadible Dimensions")
	} else {
		x1 <- subset(x,subset&!is.na(x))
		if(length(by) != 0) by1 <- subset(by,subset&!is.na(x))
	}
	if(length(x1) == 0) {
		NA
	} else {
		if(length(by) == 0) {
			func1(x1,...)
		} else {
			if(length(bys) == 0) bys <- sort(unique(by1))
			
			# there has to be a better way than this...
			dat1 <- cbind(x1,by1)[order(by1),]
			J <- length(unique(by1))
			start <- rep(NA,J)
			start[1] <- 1
			tab1 <- table(by1)
			if(J > 1) {
				for(j in 2:J) start[j] <- start[j-1] + tab1[j-1]
				end <- c(start[2:J]-1,length(x1))
			} else {
				end <- length(x)
			}
			ret <- rep(NA,length(bys))
			for(j in 1:J) {
				k <- fmatch(dat1[start[j],2],bys)
				if(!is.na(k)) ret[k] <- func1(as.numeric(dat1[start[j]:end[j],1]),...)
			}
			names(ret) <- bys
			ret
		}
	}
}

str_split_vec <- function(x,sep)
{
	temp <- str_split(x,sep)
	n <- length(temp)
	m <- length(temp[[1]])
	res <- matrix(rep(NA,n*m),ncol=m)
	for(i in 1:n)
	{
		if(length(temp[[i]]) != m) stop("error in str_split_vec -- columns are not equal size")
		for(j in 1:m) res[i,j] <- temp[[i]][j]
	}	
	res
}

mean0 <- function(x,w=numeric(0),subset=numeric(0),by=numeric(0),bys=numeric(0),se=FALSE,ci=FALSE)
{
	n <- length(x)
	if(length(subset) != 0 & length(subset) != n) stop("error in mean0 -- incompadible dimensions (subset)")
	if(length(w) != 0 & length(w) != n) stop("error in mean0 -- incompadible dimensions (w)")
	if(!se & !ci) {
		if(length(w) == 0) {
			apply0(mean,x,subset,by,bys)
		} else {
			if(length(subset) == 0) {
				subset1 <- !is.na(x) & !is.na(w)
			} else {
				subset1 <- !is.na(x) & !is.na(w) & subset
			}
			apply0(mean,w*x,subset1,by,bys) / apply0(mean,w,subset1,by,bys)
		}
	} else {
		ret <- list()
		if(length(w) == 0) {
			ret$est <- apply0(mean,x,subset,by,bys)
			se1 <- ret$se <- (var0(x,subset,by,bys) / n0(x,subset,by,bys))^.5
		} else {
			if(length(subset) == 0) {
				subset1 <- !is.na(x) & !is.na(w)
			} else {
				subset1 <- !is.na(x) & !is.na(w) & subset
			}
			W = apply0(sum,w,subset1,by,bys)
			ret <- list()
			ret$est <- apply0(sum,w*x,subset1,by,bys) / W 
			se1 <- (var0(w*x,subset1,by,bys) * n0(w*x,subset1,by,bys) / W^2)^.5
		}
		if(se) ret$se <- se1
		if(ci) ret$ci <- c(ret$est - 1.96 * se1,ret$est + 1.96 * se1)
		ret
	}
}

# sd, removing missing values, and subset option
sd0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	if(length(subset) == 0) {
		subset1 <- subset(x,!is.na(x))
		if(length(by) != 0) by1 <- subset(by,!is.na(x))
	} else if(length(subset) != length(x)) {
		# Error
		stop("Incompadible Dimensions")
	} else {
		subset1 <- subset(x,subset&!is.na(x))
		if(length(by) != 0) by1 <- subset(by,subset&!is.na(x))
	}
	if(length(subset1) == 0) {
		NA
	} else {
		if(length(by) == 0) {
			sd(subset1)
		} else {
			if(length(bys) == 0) bys <- sort(unique(by1))
			n <- length(bys)
			ret <- rep(NA,n)
			for(i in 1:n)
			{
				if(length(subset(subset1,by1==bys[i])) > 0) {
					ret[i] <- sd(subset(subset1,by1==bys[i]))
				}
			}
			names(ret) <- bys
			ret
		}
	}
}

max0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	apply0(max,x,subset,by,bys)
}

min0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	apply0(min,x,subset,by,bys)
}

sum0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	apply0(sum,x,subset,by,bys)
}

median0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0),se=FALSE)
{
	if(!se) {
		apply0(median,x,subset,by,bys)
	} else {
		list(est=apply0(median,x,subset,by,bys),se=apply0(medianse,x,subset,by,bys))
	}
}


# sample size for 1-D statistics
n0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0))
{
	if(length(subset) == 0) {
		subset1 <- subset(x,!is.na(x))
		if(length(by) != 0) by1 <- subset(by,!is.na(x))
	} else if(length(subset) != length(x)) {
		# Error
		stop("Incompadible Dimensions")
	} else {
		subset1 <- subset(x,subset&!is.na(x))
		if(length(by) != 0) by1 <- subset(by,subset&!is.na(x))
	}
	if(length(subset1) == 0) {
		0
	} else {
		if(length(by) == 0) {
			sum(!is.na(subset1))
		} else {
			if(length(bys) == 0) bys <- sort(unique(by1))
			n <- length(bys)
			ret <- rep(NA,n)
			for(i in 1:n)
			{
				if(length(subset(subset1,by1==bys[i])) > 0) {
					ret[i] <- sum(!is.na(subset(subset1,by1==bys[i])))
				}
			}
			names(ret) <- bys
			ret[is.na(ret)] <- 0
			ret
		}
	}
}

quantile0 <- function(x,probs,subset=logical(0),by=numeric(0),bys=numeric(0))
{
	n <- length(x)
	if(length(subset) != 0 & length(subset) != n) stop("error in quantile0 -- incompadible dimensions (subset)")
	k <- length(probs)
	if(length(by) == 0) {
		m <- 1
	} else {
		if(length(bys) == 0) bys <- sort(unique(by))
		m <- length(bys)
	}
	ret <- matrix(rep(NA,m*k),nrow=m)
	for(i in 1:k) ret[,i] <- apply0(quantile,x,subset,by,bys,probs=probs[i])
	rownames(ret) <- bys
	colnames(ret) <- probs
	ret
}

medianse <- function(x)
{
	m1 <- median(x)
	d1 <- density(x,from=m1,to=m1,n=1)$y
	1 / (2 * d1 * length(x)^.5)	
}

# vectorized function that determines if a row of a matrix has any missing values
rowmiss <- function(X)
{
	#!apply(X, 1, function(x)!any(is.na(x)))
	rowSums(is.na(X))>0 # faster?
}

# function for robust standard errors
robust0 <- function(fm)
{
	require(sandwich, quietly = TRUE)
	require(lmtest, quietly = TRUE)
	N <- length(fm$residuals)
	M <- N
	K <- fm$rank
	dfc <- (M/(M-1))*((N-1)/(N-K))
	uj  <- apply(estfun(fm),2, function(x) tapply(x, 1:N, sum))
	vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
	sm <- summary(fm)
	model <- coeftest(fm, vcovCL)
	list(model=model,N=N,r.squared=sm$r.squared,coef=sm$coefficients[1:K],theta=sm$coefficients[1:K],se=(diag(vcovCL))^.5,V=vcovCL)	
}

# function for clustered standard errors
cluster0 <- function(fm, cl)
{
	require(sandwich, quietly = TRUE)
	require(lmtest, quietly = TRUE)
	not <- attr(fm$model,"na.action")
	if(!is.null(not)) { cl <- cl[-not] }
	M <- length(unique(cl))
	N <- length(cl)
	K <- fm$rank
	dfc <- (M/(M-1))*((N-1)/(N-K))
	uj <- apply(estfun(fm),2, function(x) tapply(x, cl, sum));
	vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
	sm <- summary(fm)
	resid <- summary(fm)$residuals
	if(length(resid) > 0) {
		sigmasqse <- NA
		model <- coeftest(fm, vcovCL)
		list(model=model,N=N,Clusters=M,r.squared=summary(fm)$r.squared,sigmasq=summary(fm)$sigma^2,sigmasqse=sigmasqse,coef=sm$coefficients[1:K],se=diag(vcovCL)^.5,V=vcovCL)

	} else {
		model <- coeftest(fm, vcovCL)
		list(model=model,N=N,Clusters=M,r.squared=fm$r.squared,coef=fm$coefficients[1:K],se=diag(vcovCL)^.5,V=vcovCL)
	}
}

# correlation, pairwise
cor0 <- function(x,y,w=logical(0),subset=logical(0),by=logical(0),bys=logical(0),x.se=logical(0),y.se=logical(0))
{
	if(length(by) == 0)
	{
		if(length(subset) == 0) {
			if(length(w) == 0) {
				if(length(x.se) == 0) Rx <- 1
				else                  Rx <- (var0(x) - mean0(x.se^2)) / var0(x)
				if(length(y.se) == 0) Ry <- 1
				else                  Ry <- (var0(y) - mean0(y.se^2)) / var0(y)
				if(Rx <= 0) stop(paste("error in cor0 -- Rx < 0. Rx =",Rx))
				if(Rx <= 0) stop(paste("error in cor0 -- Ry < 0. Ry =",Ry))
				cor(x,y,use="pairwise") / sqrt(Rx * Ry)
			} else {
				subset <- !is.na(x) & !is.na(y) & !is.na(w)
				x0 <- subset(x,subset)
				y0 <- subset(y,subset)
				w0 <- subset(w,subset)
				if(length(x.se) > 0) x.se0 <- subset(x.se,subset)
				if(length(y.se) > 0) y.se0 <- subset(y.se,subset)
				if(length(x.se) == 0) Rx <- 1
				else                  Rx <- (var(x0) - mean0(x.se0^2)) / var(x0)
				if(length(y.se) == 0) Ry <- 1
				else                  Ry <- (var(y0) - mean0(y.se0^2)) / var(y0)
				if(Rx <= 0) stop(paste("error in cor0 -- Rx < 0. Rx =",Rx))
				if(Rx <= 0) stop(paste("error in cor0 -- Ry < 0. Ry =",Ry))
				cov.wt(cbind(x0,y0),wt=w0,cor=TRUE)$cor[1,2] / sqrt(Rx * Ry)
			}
		} else {
			if(length(w) == 0) {
				subset1 <- subset & !is.na(x) & !is.na(y)
				x0 <- subset(x,subset1)
				y0 <- subset(y,subset1)
				if(length(x.se) > 0) x.se0 <- subset(x.se,subset1)
				if(length(y.se) > 0) y.se0 <- subset(y.se,subset1)
				if(length(x.se)==0) Rx <- 1
				else                Rx <- (var(x0) - mean0(x.se0^2)) / var(x0)
				if(length(y.se)==0) Ry <- 1
				else                Ry <- (var(y0) - mean0(y.se0^2)) / var(y0)
				if(Rx <= 0) stop(paste("error in cor0 -- Rx < 0. Rx =",Rx))
				if(Rx <= 0) stop(paste("error in cor0 -- Ry < 0. Ry =",Ry))
				cor(x0,y0,use="pairwise") / sqrt(Rx * Ry)
			} else {
				subset1 <- subset & !is.na(x) & !is.na(y) & !is.na(w)
				x0 <- subset(x,subset1)
				y0 <- subset(y,subset1)
				w0 <- subset(w,subset1)
				if(length(x.se) > 0) x.se0 <- subset(x.se,subset1)
				if(length(y.se) > 0) y.se0 <- subset(y.se,subset1)
				if(length(x.se) == 0) Rx <- 1
				else                  Rx <- (var0(x0,w=w0) - mean0(x.se0^2,w=w0)) / var0(x0,w=w0)
				if(length(y.se) == 0) Ry <- 1
				else                  Ry <- (var0(y0,w=w0) - mean0(y.se0^2,w=w0)) / var0(y0,w=w0)
				if(Rx <= 0) stop(paste("error in cor0 -- Rx < 0. Rx =",Rx))
				if(Rx <= 0) stop(paste("error in cor0 -- Ry < 0. Ry =",Ry))
				cov.wt(cbind(x0,y0),wt=w0,cor=TRUE)$cor[1,2] / sqrt(Rx * Ry)
			}
		}
	} else {
		if(length(bys) == 0) {
			if(length(subset) == 0) bys <- sort(unique(by))
			else                    bys <- sort(unique(subset(by,subset)))
		}
		ret <- rep(NA,length(bys))
		for(i in 1:length(bys)) {
			if(length(subset) == 0) {
				x0 <- subset(x,by==bys[i])
				y0 <- subset(y,by==bys[i])
				if(length(x.se) > 0) x.se0 <- subset(x.se,by==bys[i])
				if(length(y.se) > 0) y.se0 <- subset(y.se,by==bys[i])
			} else {
				x0 <- subset(x,subset&by==bys[i])
				y0 <- subset(y,subset&by==bys[i])
				if(length(x.se) > 0) x.se0 <- subset(x.se,subset&by==bys[i])
				if(length(y.se) > 0) y.se0 <- subset(y.se,subset&by==bys[i])
			}
			if(length(x0) > 0) {
				
				if(length(x.se)==0) Rx <- 1
				else                Rx <- (var(x0) - mean0(x.se0^2)) / var(x0)
				if(length(y.se)==0) Ry <- 1
				else                Ry <- (var(y0) - mean0(y.se0^2)) / var(y0)
				if(Rx <= 0) stop(paste("error in cor0 -- Rx < 0. Rx =",Rx))
				if(Rx <= 0) stop(paste("error in cor0 -- Ry < 0. Ry =",Ry))
				if(length(w) == 0) {
					ret[i] <- cor(x0,y0,use="pairwise") / sqrt(Rx * Ry)
				} else {
					w0 <- subset(w,subset)
					ret[i] <- cov.wt(cbind(x0,y0),wt=w0,cor=TRUE)$cor[1,2] / sqrt(Rx * Ry)
				}
			}
		}
		names(ret) <- bys
		ret
	}
}

var0 <- function(x,subset=numeric(0),by=numeric(0),bys=numeric(0),w=numeric(0))
{
	if(length(w) == 0) {
		apply0(var,x,subset,by,bys)
	} else {
		if(length(by) > 0) {
			stop("error in var0 -- this part not done yet!")
		} else {
			# weighted variance
			if(length(subset) == 0) {
				w1 <- w[!is.na(x)] 
				x1 <- x[!is.na(x)]
			} else {
				w1 <- w[!is.na(x) & subset] 
				x1 <- x[!is.na(x) & subset]
			}
			sum.w <- sum(w1) 
			sum.w2 <- sum(w1^2) 
			mean.w <- sum(x1 * w1) / sum(w1) 
			(sum.w / (sum.w^2 - sum.w2)) * sum(w1 * (x1 - mean.w)^2)
		}
	}
}

reli0 <- function(x,x.se,subset=logical(0))
{
	(var0(x,subset=subset) - mean0(x.se^2,subset=subset)) / var0(x,subset=subset)
}

as.per <- function(x,k=1,tex=FALSE)
{
	if(tex) {
		paste(dig(100*x,k=k),"\\%",sep="")
	} else {
		paste(dig(100*x,k=k),"%",sep="")
	}
}

as.per0 <- function(x,tex=FALSE)
{
	as.per(x,k=0,tex)
}

density0 <- function(x,bw=logical(0),subset=logical(0),w=logical(0))
{
	if(length(w) == 0) {
		if(length(subset) == 0) {
			x1 <- x[!is.na(x)]
		} else {
			x1 <- x[!is.na(x)&!is.na(subset)&subset]
		}
		if(length(bw) == 0) {
			density(x1)
		} else {
			density(x1,bw)
		}
	} else {
		if(length(subset) == 0) {
			x1 <- x[!is.na(x) & !is.na(w)]
			w1 <- w[!is.na(x) & !is.na(w)]
		} else {
			x1 <- x[!is.na(x)&!is.na(w)&!is.na(subset)&subset]
			w1 <- w[!is.na(x)&!is.na(w)&!is.na(subset)&subset]
		}
		w1 <- w1 / sum(w1)
		if(length(bw) == 0) {
			density(x1,weights=w1)
		} else {
			density(x1,bw,weights=w1)
		}
	}	
}

read0 <- function(file,sep="\t")
{
	read.delim(file,stringsAsFactors=FALSE,quote="",sep=sep)
}

dig <- function(x,k=3)
{
	y <- format(round(x,k),nsmall=k)
	for(i in 1:length(y)) {
		if(is.na(x[i])) y[i] <- ""
	}
	y
}

dig1 <- function(x)
{
	dig(x,1)
}

dig2 <- function(x)
{
	dig(x,2)
}

dig3 <- function(x)
{
	dig(x,3)
}

pval <- function(est,se)
{
	2*pnorm(-abs(est/se))
}

stars.scalar <- function(coef,se)
{
	if(is.na(coef)) {
		""
	} else if(is.na(se)) {
		"?"
	} else {
		pval1 <- pval(coef,se)
		if(pval1 <= 0.001) {
			"***"
		} else if(pval1 <= 0.01) {
			"**"
		} else if(pval1 <= 0.05) {
			"*"
		} else if(pval1 <= 0.1) {
			"+"
		} else {
			""
		}
	}
}

stars <- function(coef,se)
{
	if(is.matrix(coef))
	{
		m <- nrow(coef)
		n <- ncol(coef)
		res <- coef
		for(i in 1:m) for(j in 1:n) res[i,j] <- stars.scalar(coef[i,j],se[i,j])
	} else {
		m <- length(coef)
		res <- coef
		for(i in 1:m) res[i] <- stars.scalar(coef[i],se[i])
	}
	res
}

as.coef <- function(coef,se,k=3,p.val=logical(0))
{
	if(length(p.val) == 0) {
		if(is.matrix(coef)) {
			res <- coef
			for(i in 1:ncol(coef)) res[,i] <- paste(dig(coef[,i],k),stars(coef[,i],se[,i]),sep="")
			res
		} else {
			paste(dig(coef,k),stars(coef,se),sep="")
		}
	} else {
		if(is.matrix(coef)) {
			res <- coef
			for(i in 1:ncol(coef)) res[,i] <- paste(dig(coef[,i],k),stars.pval(coef[,i],p.val[,i]),sep="")
			res
		} else {
			paste(dig(coef,k),stars.pval(p.val),sep="")
		}
	}
}

as.se <- function(se,k=3)
{
	if(is.matrix(se)) {
		res <- se
		for(i in 1:ncol(se)) res[,i] <- str_replace_all(paste("(",dig(se[,i],k),")",sep=""),"\\(\\)","")
		res
	} else {
		str_replace_all(paste("(",dig(se,k),")",sep=""),"\\(\\)","")
	}
}

# get the model matrix, including NAs, from a multipart formula
model.matrix0 <- function(f1,rhs=logical(0))
{
	mf <- model.frame(f1,na.action=na.pass)
	if(length(rhs) == 0) mm <- model.matrix(f1,data=mf)
	else                 mm <- model.matrix(f1,data=mf,rhs=rhs)
	mm
}

# get the model response, including NAs, from a multipart formula
model.response0 <- function(f1,lhs=logical(0))
{
	mf <- model.frame(f1,na.action=na.pass)
	if(length(lhs) == 0) mr <- model.response(data=mf)
	else                 mr <- model.part(f1,data=mf,lhs=lhs)
	mr
}

# get all the y's and X's from a multipart formula, eliminating NAs and subseting
formula.data <- function(f,subset=logical(0),list1=logical(0),se="",nointercept=FALSE)
{
	f1 <- Formula(f)
	len1 <- length(f1)
	if(length(len1) != 2) stop("error in formula.data -- invalid formula")
	y <- list()
	X <- list()
	for(i in 1:len1[1]) y[[i]] <- model.response0(f1,lhs=i)
	for(i in 1:len1[2]) X[[i]] <- model.matrix0(f1,rhs=i)
	ptm <- proc.time()
	y0 <- y
	X0 <- X
	n <- nrow(as.matrix(y[[1]]))
	drop <- rep(0,n)
	for(i in 1:len1[1]) drop <- drop | is.na(y[[i]])
	for(i in 1:len1[2]) drop <- drop | rowmiss(X[[i]])
	if(length(subset) == n) {
		drop <- drop | !subset
		drop[is.na(drop)] <- TRUE
	} else if(length(subset) != 0) {
		stop("error in formula data")
	}
	for(i in 1:len1[1]) y[[i]] <- y[[i]][!drop]
	for(i in 1:len1[2]) X[[i]] <- X[[i]][!drop,]
	if(nointercept) {
		names1 <- colnames(X[[i]])[2:ncol(X[[i]])]
		for(i in 1:len1[2]) X[[i]] <- as.matrix(X[[i]][,2:ncol(X[[i]])])
		colnames(X[[i]]) <- names1
	}	
	if(length(list1) > 0) for(i in 1:length(list1)) if(length(list1[[i]]) > 0)
	{
		if(is.matrix(list1[[i]])) {
			if(nrow(list1[[i]]) != n) stop("error in formula.data -- element of list1 does not have n rows")
			list1[[i]] <- as.matrix(list1[[i]][!drop,])
		} else if(is.array(list1[[i]])) {
			if(length(dim(list1[[i]])) != 3) stop("error in formula.data -- array is not a 3D array")
			if(dim(list1[[i]])[1] != n) stop("error in formula.data -- element of list1 does not have n rows")
			list1[[i]] <- as.array(list1[[i]][!drop,,])
		} else {
			if(length(list1[[i]]) != n) stop("error in formula.data -- element of list1 is not of length n")
			list1[[i]] <- list1[[i]][!drop]
		}
	}

	# code to drop clusters for panel Newey-West
	if(se == "pnw" | se == "panelneweywest")
	{
		if(length(list1$unitid) == 0) stop("error in formula.data -- panel newey west standard errors requested, but unitid missing")
		if(length(list1$timeid) == 0) stop("error in formula.data -- panel newey west standard errors requested, but timeid missing")
		N <- length(y[[1]])
		miss <- rep(FALSE,N)
		for(i in 1:len1[1]) miss <- miss | is.na(y[[i]])
		for(i in 1:len1[2]) miss <- miss | is.na(X[[i]])
		unitids <- sort(unique(list1$unitid))
		T <- max(list1$timeid) - min(list1$timeid) + 1
		lag <- floor(4 * (T / 100)^(4/9))	
		pnwuse <- rep(TRUE,N)
		for(i in 1:length(unitids))
		{
			timecurr <- list1$timeid[list1$unitid==unitids[i]]
			misscurr <- miss[list1$unitid==unitids[i]]
			for(j in 1:length(timecurr)) if(misscurr[j]) timecurr[j] <- NA
			use <- TRUE
			for(j in 1:lag) if(sum(!is.na(fmatch(timecurr-j,timecurr))) < 20) use <- FALSE
			if(!use)
			{
				print(paste("Dropping cluster: ",unitids[i],sep=""))
				pnwuse[list1$unitid==unitids[i]] <- FALSE
			}
		}
		for(i in 1:len1[1]) y[[i]] <- y[[i]][pnwuse]
		for(i in 1:len1[2]) X[[i]] <- X[[i]][pnwuse,]
		list1$timeid <- list1$timeid[pnwuse]
		list1$unitid <- list1$unitid[pnwuse]
		if(length(y[[1]]) < 20) stop("error in formula.data -- panel newey west standard errors requested, but less than 20 observations are left")
	}

	for(i in 1:len1[2]) {
		keepcol <- colSums(abs(X[[i]]))!=0 # drop columns that are all zero (typically created by factor with subset
		X[[i]] <- X[[i]][,keepcol]
		X0[[i]] <- X0[[i]][,keepcol]
	}
	list(y=y,X=X,y0=y0,X0=X0,list1=list1)
}

get.ses <- function(m1,se,unitid=logical(0),timeid=logical(0))
{
	     if(se == "classic"                     ) stop("error in get.ses --- don't call with classic")
	else if(se == "cluster"                     ) res <- cluster0(m1,unitid)
	else if(se == "nw" | se == "neweywest"      ) res <- neweywest0(m1,timeid)
	else if(se == "pnw" | se == "panelneweywest") res <- panelneweywest0(m1,timeid,unitid)
	else if(se == "robust")                       res <- robust0(m1) # the default option
	else                                          stop("invalid se option")
}

# linear regression, with subset, no intercept, and cluster options
lm0 <- function(f1,se="robust",w=logical(0),subset=logical(0),nointercept=FALSE,unitid=logical(0),timeid=logical(0))
{
	# get y and x
	dat <- formula.data(f1,subset,list(unitid=unitid,timeid=timeid,w=w),se=se,nointercept=nointercept)
	if(length(dat$y) != 1) stop("error in lm0 -- invalid number of repsonses")
	if(length(dat$X) != 1) stop("error in lm0 -- invalid RHS")
	y1 <- dat$y[[1]]
	x1 <- dat$X[[1]]
	w1 <- dat$list1$w
	unitid1 <- dat$list1$unitid
	timeid1 <- dat$list1$timeid
	
	# estimate the model
	if(length(w1) == 0) {
		lm1 <- lm(y1 ~ 0 + x1)
	} else {
		lm1 <- lm(y1 ~ 0 + x1,weights=w1)
	}
	if(se=="classic") {
		s1 <- summary(lm1)
		res <- list()
		res$coef <- s1$coefficients[,1]
		res$se <- s1$coefficients[,2]
		res$N <- length(y1)
		res$r.squared <- s1$r.squared
		res$ivs <- names(lm1$coefficients)
	} else if(se == "cluster-wild") {
		K <- lm1$rank # number of covariates
		p.val <- lm.cluster.wild(f1,logical(0),x1=x1,y1=y1,unitid=unitid1,forceNull=FALSE)
		res <- get.ses(lm1,se="cluster",unitid=unitid1)
		res$ivs <- rownames(res$model)
		res$p.val <- p.val	
	} else {
		res <- get.ses(lm1,se,timeid=timeid1,unitid=unitid1)
		res$ivs <- rownames(res$model)
	}

	if(length(unitid1) == 0) {
		res$Clusters <- NULL
	} else {
		res$Clusters <- length(unique(unitid1))
	}
	
	# model predictions
	x <- dat$X0[[1]]
	if(nointercept) x <- matrix(x[,2:ncol(x)]) # the matrix here prevent submatrix from converting to vector
	pred <- x %*% lm1$coef
	XpX <- t(x1) %*% x1
	res$pred <- pred
	res$XpX <- XpX

	# hack to get correct model fit when intercept is included as a variable
	if(length(w1) == 0) {
		lm2 <- lm(y1 ~ x1)
	} else {
		lm2 <- lm(y1 ~ x1,weights=w1)
	}

	res$theta <- res$coef
	res$r.squared <- summary(lm2)$r.squared
	res$sigma <- summary(lm2)$sigma
	res$AIC <- AIC(lm2)
	res$BIC <- BIC(lm2)
	fac <- (res$N-1)/(res$N-ncol(x))
	res$r.squared.adj <- fac * res$r.squared + 1 - fac
	res$ivs <- substr(res$ivs,3,nchar(res$ivs))
	res$ivs[res$ivs=="tercept)"] <- "Intercept"
	if(se=="classic") {
		res$V <- res$sigma^2 * solve(XpX)
		res$predse <- rep(NA,nrow(x))
		for(n in 1:nrow(x)) res$predse[n] <- res$sigma * (1 + t(x[n,]) %*% solve(XpX) %*% x[n,])^.5
		tcrit <- qt(.975,res$N-ncol(x))
		res$predlow <- res$pred - tcrit * res$predse
		res$predhigh <- res$pred + tcrit * res$predse
	}
	res
}

# model without intercept is not allowed
oprobit0 <- function(f1,se="robust",subset=logical(0),unitid=logical(0),timeid=logical(0))
{
	# get y and x
	y <- get(all.vars(f1)[1])
	x <- model.matrix0(f1)
	n <- length(y)
	k <- ncol(x)

	# get the subset of the data
	if(length(subset) == 0) {
		subset1 <- rep(TRUE,n)
	} else if(length(subset) != n) {
		stop(paste("Incompadible Dimensions, length(subset): ",length(subset),", n: ",n,sep=""))
	} else {
		subset1 <- subset
	}
	for(i in 1:n) {
		if(is.na(subset1[i])) {
			subset1[i] <- FALSE
		}
	}
	y1 <- subset(y,subset1)
	x1 <- x[(1:n)[subset1],2:k]
	if(length(unitid) == n) {
		unitid1 <- unitid[subset1]
	} else if(length(unitid) != 0) {
		stop(paste("Incompadible Dimensions, length(unitid): ",length(unitid),sep=""))
	}
	if(length(timeid) == n) {
		timeid1 <- timeid[subset1]
	} else if(length(timeid) != 0) {
		stop(paste("Incompadible Dimensions, length(timeid): ",length(timeid),sep=""))
	}

	# code to drop cluster for panel Newey-West
	if(se == "pnw" | se == "panelneweywest")
	{
		miss <- is.na(y1) | rowmiss(x1)
		unitids <- sort(unique(unitid1))
		T <- max(timeid1) - min(timeid1) + 1
		lag <- floor(4 * (T / 100)^(4/9))	
		pnwuse <- rep(TRUE,length(y1))

		for(i in 1:length(unitids))
		{
			timecurr <- timeid1[unitid1==unitids[i]]
			misscurr <- miss[unitid1==unitids[i]]
			for(j in 1:length(timecurr))
			{
				if(misscurr[j])
				{
					timecurr[j] <- NA
				}
			}
			use <- TRUE
			for(j in 1:lag)
			{
				if(sum(!is.na(fmatch(timecurr-j,timecurr))) < 20)
				{
					use <- FALSE
				}
			}
			if(!use)
			{
				print(paste("Dropping cluster: ",unitids[i],sep=""))
				pnwuse[unitid1==unitids[i]] <- FALSE
			}
		}
		
		y1 <- y1[pnwuse]
		x1 <- x1[pnwuse,]
		timeid1 <- timeid1[pnwuse]
		unitid1 <- unitid1[pnwuse]
	}

	# estimate the model
	y2 <- factor(y1)
	oprobit1 <- polr(y2 ~ x1,method = c("probit"),Hess=TRUE)
		
	# trick robust se's into working by standarding oprobit1 object
	oprobit1$coef <- c(oprobit1$coef,oprobit1$zeta)
	oprobit1$residuals <- rep(NA,oprobit1$n)
	oprobit1$rank <- length(oprobit1$coef)
	usedord <- TRUE
	
	if(se == "classic")
	{
		ret <- summary(oprobit1)
	} else if(se == "cluster") {
		ret <- cluster0(oprobit1,unitid1)
	} else if(se == "nw" | se == "newneywest") {
		ret <- neweywest0(oprobit1,timeid1)
	} else if(se == "pnw" | se == "panelneweywest") {
		ret <- panelneweywest0(oprobit1,timeid1,unitid1)
	} else {
		ret <- robust0(oprobit1) # the defauly option
	}
	ret
}

lm.me <- function(f,subset=logical(0),R=logical(0),omega=logical(0))
{
	dat <- formula.data(f,subset,list(omega=omega))
	if(length(dat$y) != 1) stop("error in lm.me -- invalid number of repsonses")
	if(length(dat$X) != 1) stop("error in lm.me -- invalid RHS")
	y <- dat$y[[1]]
	X <- dat$X[[1]]
	K <- ncol(X)	
	N <- length(y)

	if(length(R) == 0 & length(omega) == 0) stop("error in lm.me -- R and omega both missing")
	if(length(R) > 0 & length(omega) > 0) stop("error in lm.me -- R and omega both provided")
	if(length(R) == 0) R <- diag(as.matrix(var(X[,2:K]) / (var(X[,2:K]) + rowMeans(as.matrix(dat$list1$omega)))))

	# initial values
	lm1 <- lm0(y ~ X,se="robust")
	beta <- lm1$coef
	XXp <- solve(t(X) %*% X / N)
	sigmaeta2 <- c(0,diag(var(X))[2:K] * (1 - R))
	M <- solve(diag(rep(1,K)) - XXp * diag(sigmaeta2))
	coef <- as.vector(M %*% beta)
	V <- M %*% lm1$V %*% t(M)
	se <- diag(V)^.5
	r.squared <- 1 - (sum((y - X %*% beta)^2) - N * sum(sigmaeta2*coef^2)) / (sum(y^2) - N * mean(y)^2)
	model <- cbind(as.coef(coef,se),as.se(se))
	rownames(model) <- rownames(coef)
	pred <- dat$X0[[1]] %*% coef
	if(length(omega) == 0) {
		pred.se <- NA
	} else {
		sigmaeps2 <- 0
		for(n in 1:N) sigmaeps2 <- sigmaeps2 + (y[n] - t(coef) %*% X[n,])^2 / N - t(coef[2:K]) %*% diag(dat$list1$omega[n,],nrow=K-1) %*% coef[2:K] / N
		pred.se <- rep(NA,nrow(dat$X0[[1]]))
		for(n in 1:nrow(dat$X0[[1]])) pred.se[n] <- (sigmaeps2 + t(coef[2:K]) %*% diag(omega[n,],nrow=K-1) %*% coef[2:K])^.5
	}
	list(model=model,theta=coef,coef=coef,V=V,se=se,N=N,r.squared=r.squared,pred=pred,pred.se=pred.se,ivs=lm1$ivs)
}

# compute variance-covariance matrix for func1(betahat), where func1 is vector-valued
delta.method <- function(func1,beta,V,...)
{
	c_beta <- jacobian(func1,beta,method="simple",...)
	Vc <- c_beta %*% V %*% t(c_beta)
	list(est=func1(beta,...),se=sqrt(diag(Vc)),V=Vc)
}

# various functions
median.se <- function(x,x.se,subset,by,bys,R=100)
{
	m <- length(bys)
	med1 <- rep(NA,m)
	se1 <- rep(NA,m)
	temp <- rep(NA,R)
	for(i in 1:m)
	{
		x1 <- subset(x,subset&by==bys[i])
		x.se1 <- subset(x.se,subset&by==bys[i])
		med1[i] <- median0(x1)
		for(r in 1:R) temp[r] <- median0(x1+x.se1*rnorm(length(x1)))
		se1[i] <- sd0(temp)
	}
	list(est=med1,se=se1)
}

quantile.se <- function(x,x.se,prob,subset,R=100)
{
	temp <- rep(NA,R)
	x1 <- subset(x,subset)
	x.se1 <- subset(x.se,subset)
	med1 <- median0(x1)
	for(r in 1:R) temp[r] <- quantile0(x1+x.se1*rnorm(length(x1)),probs=prob)
	se1 <- sd0(temp)
	list(est=med1,se=se1)
}

state_plot <- function(Alpha,inddat)
{
	dems <- median0(Alpha,by=inddat$StateChamber,subset=inddat$Party=="D"&inddat$Chamber!="G")
	reps <- median0(Alpha,by=inddat$StateChamber,subset=inddat$Party=="R"&inddat$Chamber!="G")
	plot(0:1,type='n',xlim=c(-5,5),ylim=c(0,102))
	text(dems,1:101,col='blue',names(dems),cex=.5)
	text(reps,1:101,col='red',names(reps),cex=.5)
}

alternate <- function(x,y)
{
	if(is.matrix(x)) {
		m <- nrow(x)
		n <- ncol(x)
		res <- matrix(rep(NA,length(x)*2),ncol=n)
		for(i in 1:n) {
			res[,i] <- as.vector(t(matrix(c(x[,i],y[,i]),nrow=m)))
		}
		res
	} else {
		n <- length(x)
		if(length(y) != n) {
			error("Error with alternate --- incompadible dimensions")
		}
		as.vector(t(matrix(c(x,y),nrow=n)))
	}
}

output.models <- function(ms)
{
	J <- length(ms)
	ivs <- ms[[1]]$ivs
	if(J > 1) for(j in 2:J) ivs <- c(ivs,ms[[j]]$ivs)
	ivs <- unique(ivs)
	K <- length(ivs)
	coef <- matrix(rep(NA,J*K),ncol=J)
	se <- matrix(rep(NA,J*K),ncol=J)
	N <- rep(NA,J)
	r.squared <- rep(NA,J)
	Clusters <- rep(NA,J)
	for(j in 1:J)
	{
		for(k in 1:length(ms[[j]]$ivs))
		{
			coef[fmatch(ms[[j]]$ivs[k],ivs),j] <- ms[[j]]$theta[k]
			se[fmatch(ms[[j]]$ivs[k],ivs),j] <- ms[[j]]$se[k]
		}
		N[j] <- ms[[j]]$N
		if(!is.null(ms[[j]]$r.squared)) r.squared[j] <- ms[[j]]$r.squared
		if(!is.null(ms[[j]]$Clusters)) Clusters[j] <- ms[[j]]$Clusters
	}
	rownames(coef) <- ivs
	colnames(coef) <- names(ms)
	rownames(se) <- ivs
	colnames(se) <- names(ms)
	res <- list(coef=coef,se=se,N=N)
	if(sum(!is.na(r.squared)) > 0) res$r.squared=r.squared
	if(sum(!is.na(Clusters)) > 0) res$Clusters=Clusters
	res
}

models.to.latex <- function(ms,clusters=TRUE)
{
	os <- output.models(ms)
	J <- nrow(os$coef)
	K <- ncol(os$coef)
	hasrs <- !is.null(os$r.squared)
	hasclusters <- !is.null(os$Clusters)
	if(!clusters) hasclusters <- FALSE
	mat1 <- matrix(rep(NA,(2*J+1+hasrs+hasclusters)*K),ncol=K)
	if(hasrs & !hasclusters) {
		for(i in 1:K) mat1[,i] <- c(alternate(as.coef(os$coef[,i],os$se[,i]),as.se(os$se[,i])),os$N[i],dig3(os$r.squared[i]))
		rownames(mat1) <- c(alternate(rownames(os$coef),rep("",J)),"N","r.squared")
	} else if(!hasrs & !hasclusters) {
		for(i in 1:K) mat1[,i] <- c(alternate(as.coef(os$coef[,i],os$se[,i]),as.se(os$se[,i])),os$N[i])
		rownames(mat1) <- c(alternate(rownames(os$coef),rep("",J)),"N")
	} else if(hasrs & hasclusters) {
		for(i in 1:K) mat1[,i] <- c(alternate(as.coef(os$coef[,i],os$se[,i]),as.se(os$se[,i])),os$N[i],os$Clusters[i],dig3(os$r.squared[i]))
		rownames(mat1) <- c(alternate(rownames(os$coef),rep("",J)),"N","Clusters","r.squared")
	} else {
		for(i in 1:K) mat1[,i] <- c(alternate(as.coef(os$coef[,i],os$se[,i]),as.se(os$se[,i])),os$N[i],os$Clusters[i])
		rownames(mat1) <- c(alternate(rownames(os$coef),rep("",J)),"N","Clusters")
	}
	colnames(mat1) <- colnames(os)
	mat1
}

latex_add_line <- function(file1,s)
{
	cat(s,file=file1,sep="\n",append=TRUE)
}

latex_table_entry <- function(file1,name,coef,bold=FALSE,italic=FALSE)
{
	if(bold)
	{
		name <- paste("\\textbf{",name,"}",sep="")
		coef <- paste("\\textbf{",coef,"}",sep="")
	}
	if(italic)
	{
		name <- paste("\\textit{",name,"}",sep="")
		coef <- paste("\\textit{",coef,"}",sep="")
	}
	s <- name
	n <- length(coef)
	for(i in 1:n) if(is.na(coef[i])) { coef[i] <- "" }
	if(n > 1) for(i in 1:(n-1)) { s <- paste(s," & ",coef[i],sep="") }
	s <- paste(s," & ",coef[n]," \\\\ \n",sep="")
	s
	cat(s,file=file1,sep="\n",append=TRUE)
}

latex_table_entries <- function(file1,names,coefs,bold=FALSE,italic=FALSE)
{
	m <- length(names)
	for(j in 1:m) latex_table_entry(file1,names[j],coefs[j,],bold)
}

latex_doc_start <- function(file1,packages=logical(0))
{
	cat("\\documentclass[11pt]{article}",file=file1,sep="\n")
	s <- "\\usepackage{float,amsmath,amssymb,amsthm,fullpage,setspace,natbib,enumerate,url,multirow,lscape"
	if(length(packages) > 0) for(i in 1:length(packages)) s <- paste(s,",",packages[i],sep="")
	s <- paste(s,"}",sep="")
	cat(s,file=file1,sep="\n",append=TRUE)
	#cat("\\usepackage[dvips]{graphicx,xcolor}",file=file1,sep="\n",append=TRUE)
	cat("\\onehalfspacing",file=file1,sep="\n",append=TRUE)
	cat("\\begin{document}",file=file1,sep="\n\n",append=TRUE)
}

latex_doc_end <- function(file1)
{
	cat("\\end{document}",file=file1,sep="\n",append=TRUE)
}

latex_table_start <- function(file1,k)
{
	cat("\\begin{table}[htp]\\centering\\scriptsize",file=file1,sep="\n",append=TRUE)
	cat("\\begin{singlespace}",file=file1,sep="\n",append=TRUE)
	cat("\\begin{tabular}{ l ",paste(rep("c",k),collapse=" ")," } \\hline \\hline",file=file1,sep="\n",append=TRUE)
}

latex_table_end <- function(file1,caption="",label="")
{
	cat("\\hline \\hline",file=file1,sep="\n",append=TRUE)
	cat("\\end{tabular}",file=file1,sep="\n",append=TRUE)
	cat("\\end{singlespace}",file=file1,sep="\n",append=TRUE)
	if(caption != "") cat(paste0(c("\\caption{\\footnotesize{",caption,"}}")),file=file1,sep="\n",append=TRUE)
	if(label != "") cat(paste0(c("\\label{",label,"}")),file=file1,sep="\n",append=TRUE)
	cat("\\end{table}",file=file1,sep="\n",append=TRUE)
}


if(fullrep_ideal)
{

# complete data
Y <- ipe_read_sparse(paste(dir1,"Y.dat",sep=""))
inddat <- read.delim(paste(dir1,"LegDat.dat",sep=""),header=TRUE,stringsAsFactors=FALSE)
inddat <- inddat[fmatch(Y$Nids,inddat[,1]),]
ushparty <- rep("",Y$N)
ushparty[inddat$StateChamber=="USH"&inddat$InRC==1] <- inddat$Party[inddat$StateChamber=="USH"&inddat$InRC==1]

# split up chambers
temp <- str_split(Y$Tids,"_")
split7 <- rep(NA,Y$T)
for(i in 1:Y$T) split7[i] <- toupper(temp[[i]][1])
chambers <- ipe_split_chamber_by_vote(Y,split7,inddat=inddat)

# NPAT ideal points
npat0 <- ipe_merge_chambers(list("BRIDGE"=chambers[[2]],"NPAT"=chambers[[3]]))
npat <- ipe_split_chamber_by_ind(npat0$Y,substr(npat0$Y$Nids,1,4)=="Resp",inddat=npat0$inddat)[[1]]
ipels <- ipe_start_sparse(npat$Y,1)
ipelu <- ipe_pmle_sparse(npat$Y,1,ipe=ipels)
ipel  <- ipe_normalize_by_group(ipelu,npat$inddat$Party,list(D=-1,R=1))

# state by state roll call vote ideal points
ipes <- list()
for(i in 1:length(chambers))
{
	print0("i",i)
	if(i != 2 & i != 3 & i != 104)
	{
		ipecurrs <- ipe_start_sparse(chambers[[i]]$Y,1)
		ipecurru <- ipe_pmle_sparse(chambers[[i]]$Y,1,ipe=ipecurrs)
		ipecurr  <- ipe_normalize_by_group(ipecurru,chambers[[i]]$inddat$Party,list(D=-1,R=1))
		ipes[[i]] <- ipecurr
	}
}

# ideal points using bridges only
bridge0 <- ipe_split_chamber_by_ind(chambers[[2]]$Y,chambers[[2]]$inddat$National=="S",inddat=chambers[[2]]$inddat)[[2]]
bridge <- ipe_split_chamber_by_vote(bridge0$Y,bridge0$Y$Tids=="BRIDGE_Favor death penalty",inddat=bridge0$inddat)[[1]]
set.seed(12345)
AlphaStart <- 0.001 * matrix(rnorm(bridge$Y$N),ncol=1)
DeltaStart <- 0.001 * matrix(rnorm(bridge$Y$T*2),ncol=2)
ipebu <- ipe_pmle_sparse(bridge$Y,1,AlphaStart=AlphaStart,DeltaStart=DeltaStart)
ipeb <- ipe_normalize_by_group(ipebu,bridge$inddat$Party,list(D=-1,R=1))

YB <- ipe_sparse_to_dense(bridge$Y)
YB[YB==0] <- NA
YB <- YB - 1
YB[,2] <- 1 - YB[,2]
YB[,10] <- 1 - YB[,10]
YB[,11] <- 1 - YB[,11]
YB[,12] <- 1 - YB[,12]
YB[,13] <- 1 - YB[,13]
YB[,14] <- 1 - YB[,14]
YB[,15] <- 1 - YB[,15]

bridgenames <- rep(NA,15)
bridgenames[1] <- "Favor school vouchers"
bridgenames[2] <- "Oppose soft money ban"
bridgenames[3] <- "Favor flat tax"
bridgenames[4] <- "Favor eliminating many business regulations"
bridgenames[5] <- "Oppose welfare spending"
bridgenames[6] <- "Missile defense spending: oppose eliminating"
bridgenames[7] <- "Missile defense spending: keep the same or increase"
bridgenames[8] <- "Missile defense spending: increase"
bridgenames[9] <- "Favor cutting taxes over strengthening soc. sec."
bridgenames[10] <- "Oppose handgun licenses"
bridgenames[11] <- "Oppose right to sue HMOs"
bridgenames[12] <- "Federal education spending: eliminate"
bridgenames[13] <- "Federal education spending: eliminate or decrease"
bridgenames[14] <- "Federal education spending: oppose increasing"
bridgenames[15] <- "Oppose prescription drug coverage for seniors"

colnames(YB) <- bridgenames

X1 <- cbind(YB,bridge$inddat$Party)
colnames(X1) <- c(colnames(YB),"Party")

# national npat ideal points
nnpat <- ipe_split_chamber_by_ind(npat$Y,npat$inddat$National,inddat=npat$inddat)[[2]]
ipenls <- ipe_start_sparse(nnpat$Y,1)
ipenlu <- ipe_pmle_sparse(nnpat$Y,1,ipe=ipenls)
ipenl  <- ipe_normalize_by_group(ipenlu,nnpat$inddat$Party,list(D=-1,R=1))
ipe_means_by_group(ipenl,nnpat$inddat$Party)
ipe_density_by_group(ipenl,nnpat$inddat$Party)

# NAES ideal points
naes0 <- ipe_merge_chambers(list("BRIDGE"=chambers[[2]],"NAES"=chambers[[104]]))
naes <- ipe_split_chamber_by_ind(naes0$Y,substr(naes0$Y$Nids,1,4)=="Resp",inddat=naes0$inddat)[[2]]
ipevs <- ipe_start_sparse(naes$Y,1)
ipevu <- ipe_pmle_sparse(naes$Y,1,ipe=ipevs)
ipev  <- ipe_normalize_by_group(ipevu,naes$inddat$Party,list(d=-1,r=1))
ipe_means_by_group(ipev,naes$inddat$Party)
ipe_density_by_group(ipev,naes$inddat$Party)

# National NPAT + voters
nnpat_naes <- ipe_merge_chambers(list("NNPAT"=nnpat,"NAES"=naes))

AlphaFac_nnpat_naes <- rep(1,nnpat_naes$Y$N)
AlphaFac_nnpat_naes[substr(nnpat_naes$Y$Nids,1,4)=="Resp"] <- 1754/58363
DeltaFac_nnpat_naes <- rep(1,nnpat_naes$Y$T)
DeltaFac_nnpat_naes[substr(nnpat_naes$Y$Tids,1,4)!="NAES"&substr(nnpat_naes$Y$Tids,1,4)!="BRID"] <- 71/742

AlphaNaes    <- rep(NA,nnpat_naes$Y$N)
AlphaNaes[fmatch(naes$Y$Nids,nnpat_naes$Y$Nids)] <- ipev$Alpha
AlphaNpatNat <- rep(NA,nnpat_naes$Y$N)
AlphaNpatNat[fmatch(nnpat$Y$Nids,nnpat_naes$Y$Nids)] <- ipenl$Alpha
AlphaBridge  <- rep(NA,nnpat_naes$Y$N)
AlphaBridge[fmatch(bridge$Y$Nids,nnpat_naes$Y$Nids)] <- ipeb$Alpha

ipennns <- list()
ipennns$Alpha <- matrix(rowMeans(cbind(lm0(AlphaBridge ~ AlphaNpatNat)$pred,lm0(AlphaBridge ~ AlphaNaes)$pred),na.rm=TRUE),ncol=1)
ipennns$Delta <- ipe_init_delta_sparse(nnpat_naes$Y,ipennns$Alpha,AlphaFac=AlphaFac_nnpat_naes,DeltaFac=DeltaFac_nnpat_naes)
ipennns <- ipe_normalize_prior(ipennns)
ipe_means_by_group(ipennns,nnpat_naes$inddat$Party) # OK
ipennnu <- ipe_pmle_sparse(nnpat_naes$Y,1,ipe=ipennns,AlphaFac=AlphaFac_nnpat_naes,DeltaFac=DeltaFac_nnpat_naes,IndWeights=AlphaFac_nnpat_naes,ItemWeights=DeltaFac_nnpat_naes)
ipennn <- ipe_normalize_by_group(ipennnu,nnpat_naes$inddat$Party,list(D=-1,R=1))
ipe_means_by_group(ipennn,nnpat_naes$inddat$Party) # OK
ipe_density_by_group(ipennn,nnpat_naes$inddat$Party,groups=c("D","d","i","r","R"))

# All survey data (NPAT and NAES)
npat_naes <- ipe_merge_chambers(list("NPAT"=npat,"NAES"=naes))

AlphaFac_npat_naes <- rep(1,npat_naes$Y$N)
AlphaFac_npat_naes[substr(npat_naes$Y$Nids,1,4)=="Resp"] <- 1754/58363
DeltaFac_npat_naes <- rep(1,npat_naes$Y$T)
DeltaFac_npat_naes[substr(npat_naes$Y$Tids,1,4)!="NAES"&substr(npat_naes$Y$Tids,1,4)!="BRID"] <- 71/742

AlphaNpat <- rep(NA,npat_naes$Y$N)
AlphaNpat[fmatch(npat$Y$Nids,npat_naes$Y$Nids)] <- ipel$Alpha
AlphaNpatNaes0 <- rep(NA,npat_naes$Y$N)
AlphaNpatNaes0[fmatch(nnpat_naes$Y$Nids,npat_naes$Y$Nids)] <- ipennn$Alpha
AlphaNpatNaes0[is.na(AlphaNpatNaes0)] <- lm0(AlphaNpatNaes0 ~ AlphaNpat)$pred[is.na(AlphaNpatNaes0)]

ipenvs <- list()
ipenvs$Alpha <- matrix(AlphaNpatNaes0,ncol=1)
ipenvs$Delta <- ipe_init_delta_sparse(npat_naes$Y,ipenvs$Alpha,AlphaFac=AlphaFac_npat_naes,DeltaFac=DeltaFac_npat_naes)
ipenvs <- ipe_normalize_prior(ipenvs)
ipenvu <- ipe_pmle_sparse(npat_naes$Y,1,ipe=ipenvs,AlphaFac=AlphaFac_npat_naes,DeltaFac=DeltaFac_npat_naes,IndWeights=AlphaFac_npat_naes,ItemWeights=DeltaFac_npat_naes)
ipenv <- ipe_normalize_by_group(ipenvu,npat_naes$inddat$Party,list(D=-1,R=1))
ipe_means_by_group(ipenv,npat_naes$inddat$Party) # OK
ipe_density_by_group(ipenv,npat_naes$inddat$Party,groups=c("D","d","i","r","R"))

# main estimation
AlphaFac <- rep(1,Y$N)
AlphaFac[substr(Y$Nids,1,4)=="Resp"] <- 1754/58363
DeltaFac <- rep(1,Y$T)
DeltaFac[substr(Y$Tids,1,4)!="NAES"&substr(Y$Tids,1,4)!="BRID"] <- 71/742

AlphaNpat2 <- rep(NA,Y$N)
AlphaNpat2[fmatch(npat$Y$Nids,Y$Nids)] <- ipel$Alpha
AlphaNpatN2 <- rep(NA,Y$N)
AlphaNpatN2[fmatch(npat$Y$Nids,Y$Nids)] <- ipel$IndFit[,1]
AlphaRcv <- matrix(rep(NA,length(chambers)*Y$N),ncol=length(chambers))
for(i in 1:length(chambers))
{
	print0("i",i)
	if(i != 2 & i != 3 & i != 104)
	{
		AlphaRcv[fmatch(chambers[[i]]$Y$Nids,Y$Nids),i] <- ipes[[i]]$Alpha
	}
}

AlphaNpatRcv <- AlphaNpat2
for(i in 1:length(chambers))
{
	print0("i",i)
	if(i != 2 & i != 3 & i != 104)
	{
		lm1 <- lm0(AlphaNpat2 ~ AlphaRcv[,i])
		AlphaNpatRcv[is.na(AlphaNpatRcv)] <- lm1$pred[is.na(AlphaNpatRcv)]
	}
}
AlphaNpatNaes <- rep(NA,Y$N)
AlphaNpatNaes[fmatch(npat_naes$Y$Nids,Y$Nids)] <- ipenv$Alpha # the NPAT data

Alpha0 <- AlphaNpatRcv
Alpha0[is.na(Alpha0)] <- lm0(Alpha0 ~ AlphaNpatNaes)$pred[is.na(Alpha0)]

ipe1s <- list()
ipe1s$Alpha <- matrix(Alpha0,ncol=1)
ipe1s$Delta <- ipe_init_delta_sparse(Y,Alpha0,AlphaFac=AlphaFac,DeltaFac=DeltaFac,IndWeights=AlphaFac,ItemWeights=DeltaFac)
ipe1u <- ipe_pmle_sparse(Y,1,ipe=ipe1s,AlphaFac=AlphaFac,DeltaFac=DeltaFac,IndWeights=AlphaFac,ItemWeights=DeltaFac)
ipe1 <- ipe_normalize_by_group(ipe1u,inddat$Party,list(D=-1,R=1))
ipe_means_by_group(ipe1,inddat$Party)

ipe1a <- ipe1
ipe1a$Alpha[ipe1a$IndFit[,1]<20] <- NA
ipe1a <- ipe_normalize_by_group(ipe1,ushparty,list(D=-1,R=1))
Alpha <- ipe1a$Alpha
AlphaSe <- ipe1a$AlphaSe

AlphaNaes2    <- rep(NA,Y$N)
AlphaNaes2[fmatch(naes$Y$Nids,Y$Nids)] <- ipev$Alpha
AlphaNaes2[ipev$IndFit[,1]<20] <- NA
AlphaNaes2 <- ipe_normalize_by_group(AlphaNaes2,inddat$Party,list(d=-1,r=1))

AlphaNpat3 <- AlphaNpat2
AlphaNpat3[AlphaNpatN2<20] <- NA
AlphaNpat3 <- ipe_normalize_by_group(AlphaNpat3,inddat$Party,list(D=-1,R=1))

AlphaRcv2 <- rowMeans(AlphaRcv,na.rm=TRUE)
AlphaRcv2[is.nan(AlphaRcv2)] <- NA

hc <- Y$Tids[substr(Y$Tids,1,24)=="NPAT_spend_healthcares98"]
currs <- sort(hc[substr(hc,30,33)=="1998"|substr(hc,30,33)=="1999"])
splithc <- !is.na(fmatch(Y$Tids,currs))
chambers_hc <- ipe_split_chamber_by_vote(Y,splithc,inddat=inddat)
Yhc <- ipe_sparse_to_dense(chambers_hc[[2]]$Y)
Yhc[Yhc==0] <- NA
Yhc <- Yhc - 1
hcspend1 <- rowSums(Yhc,na.rm=TRUE)+1
hcspend <- rep(NA,Y$N)
hcspend[fmatch(rownames(Yhc),Y$Nids)] <- hcspend1 

X2 <- cbind(Alpha,AlphaSe,AlphaNaes2,AlphaNpat3,AlphaRcv2,hcspend,inddat)
colnames(X2) <- c("alpha","alphase","alpha_naes","alpha_npat","alpha_rcv","HCSpend",colnames(inddat))

statedata <- read_excel(paste(dir1,"orig.xlsx",sep=""),skip=2,sheet="state")
GovId <- statedata$GovId
AlphaGovBonica <- statedata$AlphaGovBonica
EWM00 <- statedata$EWM00
EWN00Se <- as.numeric(statedata$EWM00Var)^.5
EWM_DemIdeo00 <- as.numeric(statedata$EWM_DemIdeo00)
EWM_RepIdeo00 <- as.numeric(statedata$EWM_RepIdeo00)
HouseMajParty <- statedata$HouseMajParty
SenMajParty <- statedata$SenMajParty
ClotureHouse <- statedata$ClotureHouse
ClotureSen <- statedata$ClotureSen
BudgetMaj <- statedata$"Budget Maj"
TaxMaj <- statedata$"Tax Maj"
Override <- statedata$Override

# Check for same order:
statelist <- sort(unique(substr(inddat$StateChamber,1,2)))
statelist <- statelist[statelist!=""]
if(sum0(statedata$State != statelist) != 0) stop("ERROR -- State not in same order!!")
S <- length(statelist)
USs <- fmatch("US",statelist)
AKs <- fmatch("AK",statelist)
HIs <- fmatch("HI",statelist)

MedianHouse <- median0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="H",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianHouseSe <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="H",by=substr(inddat$StateChamber,1,2),bys=statelist)$se
MedianSen     <- median0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="S",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianSenSe   <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="S",by=substr(inddat$StateChamber,1,2),bys=statelist)$se

MedianHouseDem   <- median0(Alpha          ,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianHouseDemSe <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)$se
MedianHouseRep   <- median0(Alpha          ,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianHouseRepSe <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)$se
MedianSenDem     <- median0(Alpha          ,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianSenDemSe   <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)$se
MedianSenRep     <- median0(Alpha          ,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)
MedianSenRepSe   <- median.se(Alpha,AlphaSe,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)$se

MeanHouse <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="H",by=substr(inddat$StateChamber,1,2),bys=statelist)
MeanSen   <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="S",by=substr(inddat$StateChamber,1,2),bys=statelist)
MeanHouseDem <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)
MeanHouseRep <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="H"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)
MeanSenDem   <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="D",by=substr(inddat$StateChamber,1,2),bys=statelist)
MeanSenRep   <- mean0(Alpha,subset=inddat$InRC==1&inddat$Chamber=="S"&inddat$Party=="R",by=substr(inddat$StateChamber,1,2),bys=statelist)

MedianResp     <- median0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",by=inddat$State,bys=statelist)
MedianRespSe   <- median0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",by=inddat$State,bys=statelist,se=TRUE)$se
MeanResp     <- mean0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",by=inddat$State,bys=statelist)
MeanRespSe     <- mean0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",by=inddat$State,bys=statelist,se=TRUE)$se
NResp          <- n0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",by=inddat$State,bys=statelist)
MedianRespDem  <- median0(Alpha,subset=inddat$Party=="d",by=inddat$State,bys=statelist)
MeanRespDem  <- mean0(Alpha,subset=inddat$Party=="d",by=inddat$State,bys=statelist)
MedianRespRep  <- median0(Alpha,subset=inddat$Party=="r",by=inddat$State,bys=statelist)
MeanRespRep  <- mean0(Alpha,subset=inddat$Party=="r",by=inddat$State,bys=statelist)
MedianVoter <- median0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MedianVoterSe <- median0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,by=inddat$State,bys=statelist,se=TRUE)$se
MeanVoter <- mean0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MeanVoterSe <- mean0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,by=inddat$State,bys=statelist,se=TRUE)$se
NVoter <- n0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MedianVoterDem <- median0(Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MeanVoterDem <- mean0(Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
NVoterDem <- n0(Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MedianVoterRep <- median0(Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
MeanVoterRep <- mean0(Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)
NVoterRep <- n0(Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1,by=inddat$State,bys=statelist)

DemFrac <- mean0(1*(inddat$Party=="d"),subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r"),by=inddat$State,bys=statelist)
DemFracSe <- mean0(1*(inddat$Party=="d"),subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r"),by=inddat$State,bys=statelist,se=TRUE)$se
RepFrac <- mean0(1*(inddat$Party=="r"),subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r"),by=inddat$State,bys=statelist)
RepFracSe <- mean0(1*(inddat$Party=="r"),subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r"),by=inddat$State,bys=statelist,se=TRUE)$se

# correct for US
MedianResp[USs]     <- median0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")
MedianRespSe[USs]   <- median0(Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",se=TRUE)$se
MeanResp[USs]       <- mean0(  Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")
MeanRespSe[USs]     <- mean0(  Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r",se=TRUE)$se
NResp[USs]          <- n0(     Alpha,subset=inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")
MedianRespDem[USs]  <- median0(Alpha,subset=inddat$Party=="d")
MeanRespDem[USs]    <- mean0(  Alpha,subset=inddat$Party=="d")
MedianRespRep[USs]  <- median0(Alpha,subset=inddat$Party=="r")
MeanRespRep[USs]    <- mean0(  Alpha,subset=inddat$Party=="r")
MedianVoter[USs]    <- median0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1)
MedianVoterSe[USs]    <- median0(Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,se=TRUE)$se
MeanVoter[USs]      <- mean0(  Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1)
MeanVoterSe[USs]    <- mean0(  Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1,se=TRUE)$se
NVoter[USs]         <- n0(     Alpha,subset=(inddat$Party=="d"|inddat$Party=="i"|inddat$Party=="r")&inddat$LikelyVoter==1)
MedianVoterDem[USs] <- median0(Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1)
MeanVoterDem[USs]   <- mean0(  Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1)
NVoterDem[USs]      <- n0(     Alpha,subset=inddat$Party=="d"&inddat$LikelyVoter==1)
MedianVoterRep[USs] <- median0(Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1)
MeanVoterRep[USs]   <- mean0(  Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1)
NVoterRep[USs]      <- n0(     Alpha,subset=inddat$Party=="r"&inddat$LikelyVoter==1)

# pivots
HouseFilLow <-  rep(NA,S)
HouseFilLowSe <-  rep(NA,S)
HouseFilHigh <- rep(NA,S)
HouseFilHighSe <- rep(NA,S)
HouseBudgetHigh <- rep(NA,S)
HouseBudgetHighSe <- rep(NA,S)
HouseTaxHigh <- rep(NA,S)
HouseTaxHighSe <- rep(NA,S)
HouseOverLow <- rep(NA,S)
HouseOverLowSe <- rep(NA,S)
HouseOverHigh <- rep(NA,S)
HouseOverHighSe <- rep(NA,S)
SenFilLow <- rep(NA,S)
SenFilLowSe <- rep(NA,S)
SenFilHigh <- rep(NA,S)
SenFilHighSe <- rep(NA,S)
SenBudgetHigh <- rep(NA,S)
SenBudgetHighSe <- rep(NA,S)
SenTaxHigh <- rep(NA,S)
SenTaxHighSe <- rep(NA,S)
SenOverLow <- rep(NA,S)
SenOverLowSe <- rep(NA,S)
SenOverHigh <- rep(NA,S)
SenOverHighSe <- rep(NA,S)

for(s in 1:S)
{
	HouseFilLow[s]       <- quantile0(Alpha,probs=1-ClotureHouse[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseFilLowSe[s]     <- quantile.se(Alpha,AlphaSe,prob=1-ClotureHouse[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se
	HouseFilHigh[s]      <- quantile0(Alpha,probs=  ClotureHouse[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseFilHighSe[s]    <- quantile.se(Alpha,AlphaSe,prob=  ClotureHouse[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se
	HouseBudgetHigh[s]   <- quantile0(Alpha,probs=BudgetMaj[s]     ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseBudgetHighSe[s] <- quantile.se(Alpha,AlphaSe,prob=BudgetMaj[s]     ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se
	HouseTaxHigh[s]      <- quantile0(Alpha,probs=TaxMaj[s]        ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseTaxHighSe[s]    <- quantile.se(Alpha,AlphaSe,prob=TaxMaj[s]        ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se
	HouseOverLow[s]      <- quantile0(Alpha,probs=1-Override[s]    ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseOverLowSe[s]    <- quantile.se(Alpha,AlphaSe,prob=1-Override[s]    ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se
	HouseOverHigh[s]     <- quantile0(Alpha,probs=  Override[s]    ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))
	HouseOverHighSe[s]   <- quantile.se(Alpha,AlphaSe,prob=  Override[s]    ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"H",sep=""))$se

	SenFilLow[s]         <- quantile0(Alpha,probs=1-ClotureSen[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenFilLowSe[s]       <- quantile.se(Alpha,AlphaSe,prob=1-ClotureSen[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
	SenFilHigh[s]        <- quantile0(Alpha,probs=  ClotureSen[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenFilHighSe[s]      <- quantile.se(Alpha,AlphaSe,prob=  ClotureSen[s],subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
	SenBudgetHigh[s]     <- quantile0(Alpha,probs=BudgetMaj[s]   ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenBudgetHighSe[s]   <- quantile.se(Alpha,AlphaSe,prob=BudgetMaj[s]   ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
	SenTaxHigh[s]        <- quantile0(Alpha,probs=TaxMaj[s]      ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenTaxHighSe[s]      <- quantile.se(Alpha,AlphaSe,prob=TaxMaj[s]      ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
	SenOverLow[s]        <- quantile0(Alpha,probs=1-Override[s]  ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenOverLowSe[s]      <- quantile.se(Alpha,AlphaSe,prob=1-Override[s]  ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
	SenOverHigh[s]       <- quantile0(Alpha,probs=  Override[s]  ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))
	SenOverHighSe[s]     <- quantile.se(Alpha,AlphaSe,prob=  Override[s]  ,subset=inddat$InRC==1&inddat$StateChamber==paste(statelist[s],"S",sep=""))$se
}

MedianHouseMaj   <- rep(NA,S)
MedianSenMaj     <- rep(NA,S)
MedianHouseMajSe <- rep(NA,S)
MedianSenMajSe   <- rep(NA,S)

for(s in 1:S)
{
	if(HouseMajParty[s] == "D") {
		MedianHouseMaj[s] <- MedianHouseDem[s]
		MedianHouseMajSe[s] <- MedianHouseDemSe[s]
	} else if(HouseMajParty[s] == "R") {
		MedianHouseMaj[s] <- MedianHouseRep[s]
		MedianHouseMajSe[s] <- MedianHouseRepSe[s]
	}

	if(SenMajParty[s] == "D") {
		MedianSenMaj[s] <- MedianSenDem[s]
		MedianSenMajSe[s] <- MedianSenDemSe[s]
	} else if(SenMajParty[s] == "R") {
		MedianSenMaj[s] <- MedianSenRep[s]
		MedianSenMajSe[s] <- MedianSenRepSe[s]
	}
}

AlphaGov <- Alpha[fmatch(statedata$GovId,Y$Nids)]
AlphaGovSe <- AlphaSe[fmatch(statedata$GovId,Y$Nids)]
AlphaGovImp <- AlphaGov
lm1 <- lm0(AlphaGov ~ statedata$AlphaGovBonica,se="classic")
AlphaGovImp[is.na(AlphaGov)] <- lm1$pred[is.na(AlphaGovImp)]
AlphaGovImpSe <- AlphaGovSe
AlphaGovImpSe[is.na(AlphaGov)] <- lm1$predse[is.na(AlphaGov)]

MedianRespImp <- MedianResp
MedianRespImpSe <- MedianRespSe
lm1 <- lm.me(MedianResp ~ statedata$EWM00,subset=(1:S)!=USs,omega=matrix(as.numeric(statedata$EWM00Var),ncol=1))
MedianRespImp[is.na(MedianRespImp)] <- lm1$pred[is.na(MedianRespImp)]
MedianRespImpSe[is.na(MedianRespImp)] <- lm1$pred.se[is.na(MedianRespImp)]
MedianRespImp[NResp<20] <- lm1$pred[NResp<20] # this does nothing!!
MedianRespImpSe[NResp<20] <- lm1$pred.se[NResp<20]

MedianVoterImp <- MedianVoter
MedianVoterImpSe <- MedianVoterSe
lm1 <- lm.me(MedianVoter ~ statedata$EWM00,subset=(1:S)!=USs,omega=matrix(as.numeric(statedata$EWM00Var),ncol=1))
MedianVoterImp[is.na(MedianRespImp)] <- lm1$pred[is.na(MedianVoterImp)]
MedianVoterImpSe[is.na(MedianRespImp)] <- lm1$pred.se[is.na(MedianVoterImp)]
MedianVoterImp[NVoter<20] <- lm1$pred[NVoter<20] # this does nothing!!
MedianVoterImpSe[NVoter<20] <- lm1$pred.se[NVoter<20]

MedianVoterDemImp <- MedianVoterDem
MedianVoterDemImp[is.na(MedianVoterDemImp)] <- lm0(MedianVoterDem ~ as.numeric(statedata$EWM_DemIdeo00),subset=(1:S)!=USs)$pred[is.na(MedianVoterDemImp)]
MedianVoterDemImp[NVoterDem<20] <- lm0(MedianResp ~ statedata$EWM00,subset=(1:S)!=USs)$pred[NVoterDem<20]

MedianVoterRepImp <- MedianVoterRep
MedianVoterRepImp[is.na(MedianVoterRepImp)] <- lm0(MedianVoterRep ~ as.numeric(statedata$EWM_RepIdeo00),subset=(1:S)!=USs)$pred[is.na(MedianVoterRepImp)]
MedianVoterRepImp[NVoterRep<20] <- lm0(MedianResp ~ statedata$EWM00,subset=(1:S)!=USs)$pred[NVoterRep<20]

StateHouseRcvFits <- rep(NA,51)
StateSenateRcvFits <- rep(NA,51)
for(i in 1:length(chambers))
{
	idx <- fmatch(chambers[[i]]$name,paste(statelist,"H",sep=""))
	if(!is.na(idx)) StateHouseRcvFits[idx] <- dig3(ipes[[i]]$Fit[6])
	idx <- fmatch(chambers[[i]]$name,paste(statelist,"S",sep=""))
	if(!is.na(idx)) StateSenateRcvFits[idx] <- dig3(ipes[[i]]$Fit[6])
}

bridge1 <- ipe_bridge_analysis(chambers)
colnames(bridge1$countind) <- colnames(bridge1$count)
colnames(bridge1$countitem) <- colnames(bridge1$count)
rownames(bridge1$countind) <- rownames(bridge1$count)
rownames(bridge1$countitem) <- rownames(bridge1$count)

mat1 <- matrix(rep(NA,6*51),nrow=51)
rownames(mat1) <- statelist

# number of bridges
temp <- rowSums(cbind(bridge1$countind[3,],bridge1$countind[,3]),na.rm=TRUE)
HouseBridges <- temp[fmatch(paste(statelist,"H",sep=""),names(temp))]
SenBridges <- temp[fmatch(paste(statelist,"S",sep=""),names(temp))]

HousePerCorr <- rep(NA,51)
SenPerCorr <- rep(NA,51)

for(i in 1:length(chambers))
{
	idx <- fmatch(chambers[[i]]$name,paste(statelist,"H",sep=""))
	if(!is.na(idx)) HousePerCorr[idx] <- dig3(ipes[[i]]$Fit[6])
	idx <- fmatch(chambers[[i]]$name,paste(statelist,"S",sep=""))
	if(!is.na(idx)) SenPerCorr[idx] <- dig3(ipes[[i]]$Fit[6])
}

# npat sample sizes
tab1 <- table(npat$inddat$StateChamber)
statetab1 <- substr(names(tab1),1,2)
chamtab1 <- substr(names(tab1),3,3)

HouseNpatN <- rep(NA,51)
SenNpatN <- rep(NA,51)
GovNpatN <- rep(NA,51)

# number of bridges
for(i in 1:length(statelist))
{
	idx <- fmatch(paste(statelist[i],"H",sep=""),names(tab1))
	if(!is.na(idx)) HouseNpatN[i] <- tab1[idx]
	idx <- fmatch(paste(statelist[i],"S",sep=""),names(tab1))
	if(!is.na(idx)) SenNpatN[i] <- tab1[idx]
	idx <- fmatch(paste(statelist[i],"G",sep=""),names(tab1))
	if(!is.na(idx)) GovNpatN[i] <- tab1[idx]
}

X3 <- cbind(statelist,MedianHouse,MedianHouseSe,MedianSen,MedianSenSe,MedianHouseDem,MedianHouseRep,MedianSenDem,MedianSenRep,MedianResp,MedianRespSe,NResp,MedianRespDem,MedianRespRep,NVoter,MedianVoter,MedianVoterSe,MedianVoterDem,MedianVoterRep,NVoterDem,NVoterRep,MeanHouse,MeanSen,MeanHouseDem,MeanHouseRep,MeanSenDem,MeanSenRep,MeanResp,MeanRespSe,MeanRespDem,MeanRespRep,MeanVoter,MeanVoterSe,MeanVoterDem,MeanVoterRep,HouseFilLow,HouseFilLowSe,HouseFilHigh,HouseFilHighSe,HouseBudgetHigh,HouseBudgetHighSe,HouseTaxHigh,HouseTaxHighSe,HouseOverLow,HouseOverLowSe,HouseOverHigh,HouseOverHighSe,SenFilLow,SenFilLowSe,SenFilHigh,SenFilHighSe,SenBudgetHigh,SenBudgetHighSe,SenTaxHigh,SenTaxHighSe,SenOverLow,SenOverLowSe,SenOverHigh,SenOverHighSe,MedianHouseMaj,MedianHouseMajSe,MedianSenMaj,MedianSenMajSe,AlphaGov,AlphaGovSe,AlphaGovImp,AlphaGovImpSe,MedianRespImp,MedianRespImpSe,MedianVoterImp,MedianVoterImpSe,MedianVoterDemImp,MedianVoterRepImp,DemFrac,DemFracSe,RepFrac,RepFracSe,StateHouseRcvFits,StateSenateRcvFits,HouseBridges,SenBridges,HouseNpatN,SenNpatN,GovNpatN,HousePerCorr,SenPerCorr)
colnames(X3) <- c("State","MedianHouse","MedianHouseSe","MedianSen","MedianSenSe","MedianHouseDem","MedianHouseRep","MedianSenDem","MedianSenRep","MedianResp","MedianRespSe","NResp","MedianRespDem","MedianRespRep","NVoter","MedianVoter","MedianVoterSe","MedianVoterDem","MedianVoterRep","NVoterDem","NVoterRep","MeanHouse","MeanSen","MeanHouseDem","MeanHouseRep","MeanSenDem","MeanSenRep","MeanResp","MeanRespSe","MeanRespDem","MeanRespRep","MeanVoter","MeanVoterSe","MeanVoterDem","MeanVoterRep","HouseFilLow","HouseFilLowSe","HouseFilHigh","HouseFilHighSe","HouseBudgetHigh","HouseBudgetHighSe","HouseTaxHigh","HouseTaxHighSe","HouseOverLow","HouseOverLowSe","HouseOverHigh","HouseOverHighSe","SenFilLow","SenFilLowSe","SenFilHigh","SenFilHighSe","SenBudgetHigh","SenBudgetHighSe","SenTaxHigh","SenTaxHighSe","SenOverLow","SenOverLowSe","SenOverHigh","SenOverHighSe","MedianHouseMaj","MedianHouseMajSe","MedianSenMaj","MedianSenMajSe","AlphaGov","AlphaGovSe","AlphaGovImp","AlphaGovImpSe","MedianRespImp","MedianRespImpSe","MedianVoterImp","MedianVoterImpSe","MedianVoterDemImp","MedianVoterRepImp","DemFrac","DemFracSe","RepFrac","RepFracSe","StateHouseRcvFits","StateSenateRcvFits","HouseBridges","SenBridges","HouseNpatN","SenNpatN","GovNpatN","HousePerCorr","SenPerCorr")

# status quo estimation
sq_npat_key <- read0(paste(dir1,"sq_NPAT_key.dat",sep=""))
NpatItems <- sq_npat_key$item
Category <- sq_npat_key$category
Wording <- sq_npat_key$wording
SqFree <- sq_npat_key$sqfree
T1 <- nrow(sq_npat_key)

data3 <- read0(paste(dir1,"sq_NPAT.dat",sep=""))
N3 <- nrow(data3)
K3 <- ncol(data3)
LegId3 <- paste("Cand_",data3$caseid,sep="")
National3 <- data3$national
Chamber3 <- data3$chamber
Dist3 <- data3$dist
First3 <- data3$firstname
Last3 <- data3$surname
Party3 <- data3$partyname
State3 <- data3$state
Year3 <- data3$year
StateChamber3 <- data3$statechamber
StateYear3 <- data3$stateyear
Alpha3 <- Alpha[fmatch(LegId3,Y$Nids)]

SpendTaxSet <- sort(unique(paste(str_split_vec(NpatItems[substr(NpatItems,1,5)=="spend"|substr(NpatItems,1,3)=="tax"],sep="_")[,1],str_split_vec(NpatItems[substr(NpatItems,1,5)=="spend"|substr(NpatItems,1,3)=="tax"],sep="_")[,2],sep="_")))
SpendTaxSet2 <- str_split_vec(sort(unique(paste(str_split_vec(NpatItems[substr(NpatItems,1,5)=="spend"|substr(NpatItems,1,3)=="tax"],sep="_")[,1],str_split_vec(NpatItems[substr(NpatItems,1,5)=="spend"|substr(NpatItems,1,3)=="tax"],sep="_")[,2],Wording[substr(NpatItems,1,5)=="spend"|substr(NpatItems,1,3)=="tax"],sep="_"))),sep="_")[,3]

KSpendTax <- length(SpendTaxSet)

SpendTax <- matrix(rep(NA,N3*KSpendTax),nrow=N3)
SpendTax3 <- matrix(rep(NA,N3*KSpendTax),nrow=N3)

for(k in 1:KSpendTax)
{
	Temp1 <- data3[,fmatch(paste(SpendTaxSet[k],"_1",sep=""),names(data3))]
	Temp2 <- data3[,fmatch(paste(SpendTaxSet[k],"_2",sep=""),names(data3))]
	Temp3 <- data3[,fmatch(paste(SpendTaxSet[k],"_3",sep=""),names(data3))]
	Temp4 <- data3[,fmatch(paste(SpendTaxSet[k],"_4",sep=""),names(data3))]
	Temp5 <- data3[,fmatch(paste(SpendTaxSet[k],"_5",sep=""),names(data3))]

	# 66666 = Greatly increase, 16666 = Slightly increase, 11666 - Maintain status, 11166 = Slightly decrease, 11116 = Greatly decrease, 11111 = Eliminate

	for(n in 1:N3)
	{
		if(Temp1[n] != 9)
		{
			if(Temp1[n] == 1 && Temp2[n] == 1 && Temp3[n] == 1) {
				SpendTax[n,k] <- 2 # Slightly decrease, Greatly decrease, Eliminate
			} else if(Temp1[n] == 1 && Temp2[n] == 1) {
				SpendTax[n,k] <- 1 # Maintain status
			} else {
				SpendTax[n,k] <- 0 # Greatly increase, Slightly increase
			}

			if(Temp5[n] == 1) {
				SpendTax3[n,k] <- 6 # Eliminate
			} else if(Temp4[n] == 1) {
				SpendTax3[n,k] <- 5 # Greatly decrease
			} else if(Temp3[n] == 1) {
				SpendTax3[n,k] <- 4 # Slightly decrease
			} else if(Temp2[n] == 1) {
				SpendTax3[n,k] <- 3 # Maintain status
			} else if(Temp1[n] == 1) {
				SpendTax3[n,k] <- 2 # Slightly increase
			} else {
				SpendTax3[n,k] <- 1 # Greatly increase
			}
		}
	}
}

stateyears <- sort(unique(StateYear3[as.integer(substr(StateYear3,3,6))>=1998]))

func1 <- function(x) { (x[2] + x[3]) / (2 * x[1]) }

curr <- 1
SqKey <- list()
SqName <- list()
SqN <- list()
Beta <- list()
BetaSe <- list()
Sq <- list()
SqSe <- list()
for(k in 1:KSpendTax)
{
	print0("k",k)
	print0("SpendTaxSet[k]",SpendTaxSet[k])

	for(i in 1:length(stateyears))
	{
		DepVar  <- SpendTax[,k][StateYear3==stateyears[i]]
		DepVar2 <- SpendTax3[,k][StateYear3==stateyears[i]]
		IndVar  <- Alpha3[StateYear3==stateyears[i]]

		if(n0(DepVar & IndVar) > 15)
		{
			if(min0(DepVar) == 0 && max0(DepVar) == 2 && length(table(DepVar)) == 3 && min(table(DepVar)) >= 2)
			{
				t <- try(o1 <- oprobit0(DepVar ~ IndVar))
				if(is.list(t))
				{
					if(!is.null(t$coef))
					{
						SqKey[[curr]] <- paste(SpendTaxSet[k],"_",stateyears[i],sep="")
						SqName[[curr]] <- paste(SpendTaxSet2[k],"_",stateyears[i],sep="")
						SqN[[curr]] <- n0(DepVar)
						dm1 <- delta.method(func1,o1$coef,o1$V)
						Beta[[curr]] <- o1$coef[1]
						BetaSe[[curr]] <- o1$se[1]
						Sq[[curr]] <- dm1$est
						SqSe[[curr]] <- dm1$se
						curr <- curr + 1
					}
				}
			}
		}
	}
}

temp <- str_split_vec(SqKey,sep="_")
temp <- temp[as.integer(substr(temp[,3],3,6))>=1998,]
Obs <- sort(unique(paste(temp[,1],temp[,2],substr(temp[,3],1,2),sep="_")))
ObsN <- length(Obs)

ObsSq <- rep(NA,ObsN)
ObsSqSe <- rep(NA,ObsN)
ObsSqBeta <- rep(NA,ObsN)
ObsSqBetaSe <- rep(NA,ObsN)
ObsSqN <- rep(NA,ObsN)
ObsOutcome <- rep(NA,ObsN)
ObsOutcomeSe <- rep(NA,ObsN)
ObsOutcomeBeta <- rep(NA,ObsN)
ObsOutcomeBetaSe <- rep(NA,ObsN)
ObsOutcomeN <- rep(NA,ObsN)
ObsState <- rep(NA,ObsN)
ObsMedVoter <- rep(NA,ObsN)

SqYear <- as.integer(substr(str_split_vec(SqKey,sep="_")[,3],3,6))
SqState <- substr(str_split_vec(SqKey,sep="_")[,3],1,2)

j <- fmatch(paste(str_split_vec(SqKey,sep="_")[,1],str_split_vec(SqKey,sep="_")[,2],substr(str_split_vec(SqKey,sep="_")[,3],1,2),sep="_"),Obs)
for(i in 1:length(SqKey))
{
	print0("i",i)
	if(SqYear[i] >= 1998 && SqYear[i] <= 2003) 
	{
		if(SqYear[i] == 1998 || SqYear[i] == 1999)
		{
			ObsSq[j[i]] = Sq[[i]];
			ObsSqSe[j[i]] = SqSe[[i]];
			ObsSqBeta[j[i]] = Beta[[i]];
			ObsSqBetaSe[j[i]] = BetaSe[[i]];
			ObsSqN[j[i]] = SqN[[i]];
		}
		else
		{
			ObsOutcome[j[i]] = Sq[[i]];
			ObsOutcomeSe[j[i]] = SqSe[[i]];
			ObsOutcomeBeta[j[i]] = Beta[[i]];
			ObsOutcomeBetaSe[j[i]] = BetaSe[[i]];
			ObsOutcomeN[j[i]] = SqN[[i]];
		}
		ObsState[j[i]] <- SqState[i]
	}
}

ObsMedVoter <- MedianVoter[fmatch(ObsState,statelist)]

X4 <- cbind(Obs,ObsState,ObsMedVoter,ObsSq,ObsSqSe,ObsSqBeta,ObsSqBetaSe,ObsSqN,ObsOutcome,ObsOutcomeSe,ObsOutcomeBeta,ObsOutcomeBetaSe,ObsOutcomeN)
colnames(X4) <- c("key","state","medvoter","sq","sq_se","sq_beta","sq_beta_se","sq_N","x","x_se","x_beta","x_beta_se","x_N")

wb <- createWorkbook("cs")
addWorksheet(wb, "bridge votes")
writeData(wb, sheet = 1,X1)
addWorksheet(wb, "ideal points")
writeData(wb, sheet = 2,X2)
addWorksheet(wb, "state data")
writeData(wb, sheet = 3,X3)
addWorksheet(wb, "sq data")
writeData(wb, sheet = 4,X4)
saveWorkbook(wb,paste(dir1,"cs.xlsx",sep=""), overwrite = TRUE)

} # end ideal

# load the data #
File1 <- read_excel(paste(dir1,"cs.xlsx",sep=""),sheet="sq data")
File2 <- read_excel(paste(dir1,"cs.xlsx",sep=""),sheet="state data")
File3 <- read_excel(paste(dir1,"cs.xlsx",sep=""),sheet="bridge votes")
File4 <- read_excel(paste(dir1,"cs.xlsx",sep=""),sheet="ideal points",col_types=c("numeric","numeric","numeric","numeric","numeric","numeric","text","text","text","text","text","text","text","numeric","text","text","numeric","text","numeric"))
File5 <- read_excel(paste(dir1,"orig.xlsx",sep=""),sheet="state_year")
File6 <- read_excel(paste(dir1,"orig.xlsx",sep=""),sheet="state",skip=2)
File7 <- read_excel(paste(dir1,"orig.xlsx",sep=""),sheet="ideal")

# state data
k2 <- nrow(File2)
state2      <- File2$State
statelong <- File6$StateLong
medvoter2   <- as.numeric(File2$MedianVoterImp)
medvoterse2 <- as.numeric(File2$MedianVoterImpSe)
gov2   <- as.numeric(File2$AlphaGovImp)
govse2 <- as.numeric(File2$AlphaGovImpSe)
medlegH2 <- as.numeric(File2$MedianHouse)
medlegS2      <- as.numeric(File2$MedianSen)
medlegHse2 <- as.numeric(File2$MedianHouseSe)
medlegSse2 <- as.numeric(File2$MedianSenSe)
rvs00 <- File6$RVS00
rvs04 <- File6$RVS04
ewm98 <- File6$EWM98
ewm99 <- File6$EWM99
ewm00 <- File6$EWM00
ewm01 <- File6$EWM01
ewm02 <- File6$EWM02
ewm03 <- File6$EWM03
ewm98se <- as.numeric(File6$EWM98Var)^.5
ewm99se <- as.numeric(File6$EWM99Var)^.5
ewm00se <- as.numeric(File6$EWM00Var)^.5
ewm01se <- as.numeric(File6$EWM01Var)^.5
ewm02se <- as.numeric(File6$EWM02Var)^.5
ewm03se <- as.numeric(File6$EWM03Var)^.5
naes_libcon <- as.numeric(File6$NAESId00)
naes_fac1 <- as.numeric(File6$NAESFac00)
cw_mass_soc  <- as.numeric(File6$cw_mass_soc)
cw_mass_econ <- as.numeric(File6$cw_mass_econ)
cw_pol_soc   <- as.numeric(File6$cw_pol_soc)
cw_pol_econ  <- as.numeric(File6$cw_pol_econ)
grum_tax <- as.numeric(File6$grum_tax)
HouseBridges <- File2$HouseBridges
SenBridges <- File2$SenBridges
HousePerCorr <- File2$HousePerCorr
SenPerCorr <- File2$SenPerCorr
HouseNpatN <- File2$HouseNpatN
SenNpatN <- File2$SenNpatN
GovNpatN <- File2$GovNpatN
st_ideo <- rowMeans(cbind(medlegH2,medlegS2,gov2))
south2 <- state2=="FL"|state2=="MS"|state2=="AL"|state2=="LA"|state2=="TX"|state2=="TN"|state2=="VA"|state2=="AR"|state2=="SC"|state2=="NC"|state2=="WV"|state2=="OK"|state2=="KY"|state2=="MD"|state2=="GA"

# policy data
k <- length(File1$state) # number of sqs
key <- File1$key
state <- File1$state
notus <- state!="US"
states_order <- c(sort(unique(state[state!="US"])),"US")
stcode <- rep(NA,k)
for(i in 1:k) stcode[i] <- which(states_order==state[i])

x         <- as.numeric(File1$x) # policy outcome
x_se      <- as.numeric(File1$x_se) # policy outcome standard error
s         <- as.numeric(File1$sq) # status quo
s_se      <- as.numeric(File1$sq_se) # status quo standard error
x_beta    <- as.numeric(File1$x_beta)
x_beta_se <- as.numeric(File1$x_beta_se)
s_beta    <- as.numeric(File1$sq_beta)
s_beta_se <- as.numeric(File1$sq_beta_se)

# drop estimates with no standard errors and with insignificant betas
for(i in 1:k)
{
	if(is.na(x_beta_se[i]) | is.na(x_se[i])) {
		x[i] <- NA
	} else if(abs(x_beta[i] / x_beta_se[i]) <= 1.645) {
		x[i] <- NA
	}
	if(is.na(s_beta_se[i]) | is.na(s_se[i])) {
		s[i] <- NA
	} else if(abs(s_beta[i] / s_beta_se[i]) <= 1.645) {
		s[i] <- NA
	}
}

# get rid of U.S.
for(i in 1:k)
{
	if(state[i] == "US")
	{
		x[i] <- NA
		s[i] <- NA
		x_se[i] <- NA
		s_se[i] <- NA
	}
}

issue <- substr(File1$key,1,nchar(File1$key)-5) # issue name
tax <- 1*(substr(issue,1,1)=="t") # tax issue dummy
issues <- sort(unique(subset(issue,state!="US")))
m <- length(issues)

medvoter <- medvoter2[fmatch(state,state2)]
medvoterse <- medvoterse2[fmatch(state,state2)]
usmed <- medvoter[state=="US"][1] # us overall median

# policy index for each state
x_index <- rep(k2,0)
x_tax_index <- rep(k2,0)
x_spend_index <- rep(k2,0)
x_index_se <- rep(k2,0)
x_tax_index_se <- rep(k2,0)
x_spend_index_se <- rep(k2,0)
for(i in 1:k2)
{
	x_index[i]           <- mean0(x,w=1/x_se^2,subset=state==state2[i])
	x_index_se[i]       <- (sum0(1 / x_se^2,subset=state==state2[i]))^-.5
	x_tax_index[i]       <- mean0(x,w=1/x_se^2,subset=state==state2[i]&tax==1)
	x_tax_index_se[i]   <- (sum0(1 / x_se^2,subset=state==state2[i]&tax==1))^-.5
	x_spend_index[i]     <- mean0(x,w=1/x_se^2,subset=state==state2[i]&tax==0)
	x_spend_index_se[i] <- (sum0(1 / x_se^2,subset=state==state2[i]&tax==0))^-.5
}

# individual level data
alpha4 <- File4$alpha
alphase4 <- File4$alphase
level4 <- File4$National
state4 <- File4$State
state4[level4=="N"] <- "US" # set state to US for national politicians
chamber4 <- File4$Chamber
party4 <- File4$Party
id4 <- File4$LegId
alpha_npat <- File4$alpha_npat
alpha_rcv <- File4$alpha_rcv
id7 <- File7$LegId
sm7 <- as.numeric(File7$sm_pos)
sm4 <- sm7[fmatch(id4,id7)]
bpr7 <- as.numeric(File7$bpr_pos)
bpr4 <- bpr7[fmatch(id4,id7)]

# state-year data
cont5 <- File5$cont
state5 <- File5$stabr
year5 <- File5$year

# control of state government
cont2 <- -as.numeric(cont5[fmatch(paste(state2,"_2000",sep=""),paste(state5,year5,sep="_"))])
temp1 <- matrix(rep(NA,20*51),ncol=20)
for(y in 1981:2000) temp1[,y-1980] <- -as.numeric(cont5[fmatch(paste(state2,y,sep="_"),paste(state5,year5,sep="_"))])
cont2 <- -as.numeric(cont5[fmatch(paste(state2,"_2000",sep=""),paste(state5,year5,sep="_"))])
contavg2 <- rowMeans(temp1)

# correlation between npat and roll call ideal points
HouseRho <- rep(NA,51)
SenRho <- rep(NA,51)
for(i in 1:51)
{
	if(state2[i] != "NE") HouseRho[i] <- dig3(cor0(alpha_npat,alpha_rcv,subset=state4==state2[i]&chamber4=="H"))
	SenRho[i] <- dig3(cor0(alpha_npat,alpha_rcv,subset=state4==state2[i]&chamber4=="S"))
}

# simple libcon scale
YB <- matrix(as.numeric(as.matrix(File3)[,1:15]),ncol=15)
libcon <- rowMeans(YB,na.rm=TRUE)

# finally done with setup, create the figures and tables for the paper

### FIGURE 1 ###
spendhealth <- File4$HCSpend
year4 <- File4$Year
sqplot <- function(stateA,stateB,type=1)
{
	stateA2 <- statelong[fmatch(stateA,state2)]
	stateB2 <- statelong[fmatch(stateB,state2)]

	alpha_A <- subset(alpha4,state4==stateA&level4=="S"&(year4==1998|year4==1999))
	alpha_B <- subset(alpha4,state4==stateB&level4=="S"&(year4==1998|year4==1999))

	spendhealth_A <- subset(spendhealth,state4==stateA&level4=="S"&(year4==1998|year4==1999))
	spendhealth_B <- subset(spendhealth,state4==stateB&level4=="S"&(year4==1998|year4==1999))

	alpha_A_1 <- mean0(spendhealth==1,subset=state4==stateA)
	alpha_A_2 <- mean0(spendhealth==2,subset=state4==stateA)
	alpha_A_3 <- mean0(spendhealth==3,subset=state4==stateA)
	alpha_A_4 <- mean0(spendhealth==4,subset=state4==stateA)
	alpha_A_5 <- mean0(spendhealth==5,subset=state4==stateA)
	alpha_A_6 <- mean0(spendhealth==6,subset=state4==stateA)

	alpha_B_1 <- mean0(spendhealth==1,subset=state4==stateB)
	alpha_B_2 <- mean0(spendhealth==2,subset=state4==stateB)
	alpha_B_3 <- mean0(spendhealth==3,subset=state4==stateB)
	alpha_B_4 <- mean0(spendhealth==4,subset=state4==stateB)
	alpha_B_5 <- mean0(spendhealth==5,subset=state4==stateB)
	alpha_B_6 <- mean0(spendhealth==6,subset=state4==stateB)

	sqA <- subset(s,key==paste("spend_healthcares98_",stateA,sep=""))
	sqB <- subset(s,key==paste("spend_healthcares98_",stateB,sep=""))

	if(type==1) {
		par(mfrow=c(2,1))
		par(mar=c(4,7, 2, 2),mgp = c(2.5, 1, 0))
		plot(alpha_A,spendhealth_A,xlim=c(-3,3),yaxt="n",ylim=c(.8,6.2),main=paste(stateA2," (s=",dig2(sqA),")",sep=""),ylab="",xlab="",col='black',font.main=1)
		axis(side=2,las=1,at=1:6,cex.axis=.8,labels=c("Lg. Increase","Sm. Increase","Maintain SQ","Sm. Decrease","Lg. Decrease","Eliminate"))
		abline(h=3,lty=1,col='grey')
		alpha_A1 <- subset(alpha_A,!is.na(spendhealth_A))
		spendhealth_A1 <- subset(spendhealth_A,!is.na(spendhealth_A))
		alpha_B1 <- subset(alpha_B,!is.na(spendhealth_B))
		spendhealth_B1 <- subset(spendhealth_B,!is.na(spendhealth_B))
		ksA <- ksmooth(alpha_A1,spendhealth_A1,bandwidth=3,n.points=512)
		ksB <- ksmooth(alpha_B1,spendhealth_B1,bandwidth=3,n.points=512)
		points(ksA$x,ksA$y,type='l',col='magenta')	
		points(ksB$x,ksB$y,type='l',col='cyan',lty=2)	

		plot(alpha_B,spendhealth_B,xlim=c(-3,3),yaxt="n",ylim=c(.8,6.2),main=paste(stateB2," (s=",round(sqB, digits=2),")",sep=""),ylab="",xlab="Ideal Point",col='black',font.main=1)
		axis(side=2,las=1,at=1:6,cex.axis=.8,labels=c("Lg. Increase","Sm. Increase","Maintain SQ","Sm. Decrease","Lg. Decrease","Eliminate"))
		abline(h=3,lty=1,col='grey')
		points(ksB$x,ksB$y,type='l',col='cyan')
		points(ksA$x,ksA$y,type='l',col='magenta',lty=2)
	} else {
		spendhealth_A2 <- spendhealth_A
		spendhealth_A2[spendhealth_A2==1|spendhealth_A2==2] <- 1
		spendhealth_A2[spendhealth_A2==3] <- 2
		spendhealth_A2[spendhealth_A2==4|spendhealth_A2==5|spendhealth_A2==6] <- 3

		spendhealth_B2 <- spendhealth_B
		spendhealth_B2[spendhealth_B2==1|spendhealth_B2==2] <- 1
		spendhealth_B2[spendhealth_B2==3] <- 2
		spendhealth_B2[spendhealth_B2==4|spendhealth_B2==5|spendhealth_B2==6] <- 3
		
		par(mfrow=c(2,1))
		par(mar=c(4,7, 2, 2),mgp = c(2.5, 1, 0))
		plot(alpha_A,spendhealth_A2,xlim=c(-3,3),yaxt="n",ylim=c(.8,3.2),main=paste(stateA2," (s=",dig2(sqA),")",sep=""),ylab="",xlab="",col='black',font.main=1)
		axis(side=2,las=1,at=1:3,cex.axis=.8,labels=c("Increase","Maintain SQ","Decrease"))
		abline(h=2,lty=1,col='grey')
		alpha_A1 <- subset(alpha_A,!is.na(spendhealth_A2))
		spendhealth_A3 <- subset(spendhealth_A2,!is.na(spendhealth_A2))
		alpha_B1 <- subset(alpha_B,!is.na(spendhealth_B2))
		spendhealth_B3 <- subset(spendhealth_B2,!is.na(spendhealth_B2))
		ksA <- ksmooth(alpha_A1,spendhealth_A3,bandwidth=3,n.points=512)
		ksB <- ksmooth(alpha_B1,spendhealth_B3,bandwidth=3,n.points=512)
		points(ksA$x,ksA$y,type='l',col='magenta')
		points(ksB$x,ksB$y,type='l',col='cyan',lty=2)

		plot(alpha_B,spendhealth_B2,xlim=c(-3,3),yaxt="n",ylim=c(.8,3.2),main=paste(stateB2," (s=",round(sqB, digits=2),")",sep=""),ylab="",xlab="Ideal Point",col='black',font.main=1)
		axis(side=2,las=1,at=1:3,cex.axis=.8,labels=c("Increase","Maintain SQ","Decrease"))
		abline(h=2,lty=1,col='grey')
		points(ksB$x,ksB$y,type='l',col='cyan')
		points(ksA$x,ksA$y,type='l',col='magenta',lty=2)
	}
}

pdf(file=paste(dir1,"fig1.pdf",sep=""),width=6,height=6)
sqplot("HI","MS",type=2)
dev.off()

### FIGURE 2 ###
pdf(file=paste(dir1,"fig2.pdf",sep=""),width=6,height=4)
plot(density0(alpha4,subset=party4=="d"),xlim=c(-3,3),ylim=c(0,1.5),col='blue',lty=2,xlab="Ideal Point",ylab="",main="")
points(density0(alpha4,subset=party4=="i"),col='black',type='l',lty=2)
points(density0(alpha4,subset=party4=="r"),col='red',type='l',lty=2)
points(density0(alpha4,subset=party4=="D"),col='blue',type='l')
points(density0(alpha4,subset=party4=="R"),col='red',type='l')
legend(-3,1.5,lty=c(1,1,2,2),col=c('blue','red','blue','black','red'),legend=c("Dem. Candidates","Rep. Candidates","Dem. Voters","Ind. Voters","Rep. Voters"),cex=.75)
dev.off()

### FIGURE 3 ###
fig3 <- matrix(rep(NA,15*5),ncol=5)
for(i in 1:15)
{
	fig3[i,1] <- mean0(YB[,i],subset=File3$Party=="D")
	fig3[i,2] <- mean0(YB[,i],subset=File3$Party=="d")
	fig3[i,3] <- mean0(YB[,i],subset=File3$Party=="d"|File3$Party=="i"|File3$Party=="r")
	fig3[i,4] <- mean0(YB[,i],subset=File3$Party=="r")
	fig3[i,5] <- mean0(YB[,i],subset=File3$Party=="R")
}

pdf(file=paste(dir1,"fig3.pdf",sep=""),width=6,height=4)
plot(0:1,type='n',xlim=c(-.6,1),ylim=c(-15,-1),ylab="",xlab="Percent Holding Conservative Position",xaxt="n",yaxt="n",cex=.75)
text(fig3[,1],-(1:15),"D",col='blue',cex=.5)
text(fig3[,2],-(1:15),"d",col='blue',cex=.5)
text(fig3[,3],-(1:15),"v",col='black',cex=.5)
text(fig3[,4],-(1:15),"r",col='red',cex=.5)
text(fig3[,5],-(1:15),"R",col='red',cex=.5)
text(-.3,-(1:15),colnames(File3)[1:15],col='black',cex=.4)
axis(1,a=c(0,.5,1), labels=c(0,.5,1))
dev.off()

### FIGURE 4 ###
pdf(file=paste(dir1,"fig4.pdf",sep=""),width=6,height=4)
mc1 <- as.matrix(read_excel(paste(dir1,"mc.xlsx",sep=""),sheet="mc1"))
cols <- c('blue','red','magenta','cyan')
par(mfrow=c(1,2))

# bias
plot(0:1,type='n',xlim=c(-2,2),ylim=c(-1,1),xlab="Alpha",ylab="Bias")
for(i in 1:4)
{
	points(lowess(mc1[1:1000,3*i-1] ~ mc1[1:1000,3*i-2]),type='l',col=cols[i])
	points(lowess(mc1[1001:1100,3*i-1] ~ mc1[1001:1100,3*i-2]),type='l',lty=2,col=cols[i])
}
legend(-2,1,c("B=5","B=10","B=20","B=100"),lty=c(1,1,1,1),col=cols,cex=.75)

# rmse
plot(0:1,type='n',xlim=c(-2,2),ylim=c(0,.4),xlab="Alpha",ylab="RMSE")
for(i in 1:4)
{
	points(lowess(mc1[1:1000,3*i] ~ mc1[1:1000,3*i-2]),type='l',col=cols[i])
	points(lowess(mc1[1001:1100,3*i] ~ mc1[1001:1100,3*i-2]),type='l',lty=2,col=cols[i])
}
dev.off()

### FIGURE 5 ###
pdf(file=paste(dir1,"fig5.pdf",sep=""),width=6,height=4)
mc2 <- as.matrix(read_excel(paste(dir1,"mc.xlsx",sep=""),sheet="mc2"))
cols <- c('blue','red','magenta','cyan')
par(mfrow=c(1,2))

# bias
plot(0:1,type='n',xlim=c(-2,2),ylim=c(-1,1),xlab="Alpha",ylab="Bias")
for(i in 1:4)
{
	points(lowess(mc2[1:10000,3*i-1] ~ mc2[1:10000,3*i-2]),type='l',col=cols[i])
	points(lowess(mc2[10001:11000,3*i-1] ~ mc2[10001:11000,3*i-2]),type='l',lty=2,col=cols[i])
}
legend(-2,1,c("B=5","B=10","B=20","B=100"),lty=c(1,1,1,1),col=cols,cex=.75)

# rmse
plot(0:1,type='n',xlim=c(-2,2),ylim=c(0,.6),xlab="Alpha",ylab="RMSE")
for(i in 1:4)
{
	points(lowess(mc2[1:10000,3*i] ~ mc2[1:10000,3*i-2]),type='l',col=cols[i])
	points(lowess(mc2[10001:11000,3*i] ~ mc2[10001:11000,3*i-2]),type='l',lty=2,col=cols[i])
}
dev.off()

### FIGURE 6 ###
pdf(file=paste(dir1,"fig6.pdf",sep=""),width=6,height=4)
mc3 <- read_excel(paste(dir1,"mc.xlsx",sep=""),sheet="mc3")

par(mfrow=c(1,3))

plot(1:5,mc3$bias[1:5],type='l',col='blue',ylim=c(-1,1),xlab="N",ylab="Bias",xaxt="n")
axis(1, at=1:5, labels=c(10,20,40,80,160))
points(1:5,mc3$bias[6:10],type='l',col='red')
points(1:5,mc3$bias[11:15],type='l',col='magenta')
points(1:5,mc3$bias[16:20],type='l',col='cyan')
points(1:5,mc3$bias[21:25],type='l',col='purple')
legend(2,1,c("sq=-2","sq=-1","sq=0","sq=1","sq=2"),lty=c(1,1,1,1,1),col=c('blue','red','magenta','cyan','purple'))

plot(1:5,mc3$rmse[1:5],type='l',col='blue',ylim=c(0,1.2),xlab="N",ylab="RMSE",xaxt="n")
axis(1, at=1:5, labels=c(10,20,40,80,160))
points(1:5,mc3$rmse[6:10],type='l',col='red')
points(1:5,mc3$rmse[11:15],type='l',col='magenta')
points(1:5,mc3$rmse[16:20],type='l',col='cyan')
points(1:5,mc3$rmse[21:25],type='l',col='purple')

plot(1:5,mc3$cov[1:5],type='l',col='blue',ylim=c(0,2),xlab="N",ylab="95% Coverage",xaxt="n")
axis(1, at=1:5, labels=c(10,20,40,80,160))
points(1:5,mc3$cov[6:10],type='l',col='red')
points(1:5,mc3$cov[11:15],type='l',col='magenta')
points(1:5,mc3$cov[16:20],type='l',col='cyan')
points(1:5,mc3$cov[21:25],type='l',col='purple')
abline(h=.95,lty=2)
dev.off()

### FIGURE 7 ### --- not based on calculations

file1 <- paste(dir1,"tables.tex",sep="")
latex_doc_start(file1)

### TABLE 1 ### --- not based on calculations

### TABLE 2 ###
tab2 <- matrix(rep(NA,8*2),ncol=2)
tab2[1,1] <- dig2(cor0(medvoter2,rvs00))
tab2[1,2] <- dig2(cor0(medvoter2,rvs00,x.se=medvoterse2))
tab2[2,1] <- dig2(cor0(medvoter2,ewm00))
tab2[2,2] <- dig2(cor0(medvoter2,ewm00,x.se=medvoterse2,y.se=ewm00se))

tab2[3,1] <- dig2(cor0(medvoter2,naes_libcon))
tab2[4,1] <- dig2(cor0(medvoter2,naes_fac1))
tab2[5,1] <- dig2(cor0(bpr4,alpha4,subset=chamber4=="H"))
tab2[6,1] <- dig2(cor0(bpr4,alpha4,subset=chamber4=="S"))
tab2[7,1] <- dig2(cor0(sm4,alpha4,subset=chamber4=="H"))
tab2[8,1] <- dig2(cor0(sm4,alpha4,subset=chamber4=="S"))

rownames(tab2) <- c("\\hspace{.4cm}Right Wing Vote Share","\\hspace{.4cm}State Liberalism \\citep{EriksonEtAl1993}","\\hspace{.4cm}NAES Self-identification","\\hspace{.4cm}NAES Principal Components \\citep{Peress13b}","\\hspace{.4cm}Lower Chambers \\citep{BattistaEtAl13}","\\hspace{.4cm}Upper Chambers \\citep{BattistaEtAl13}","\\hspace{.4cm}Lower Chambers \\citep{ShorMcCarty11}","\\hspace{.4cm}Upper Chambers \\citep{ShorMcCarty11}")

latex_table_start(file1,2)
latex_add_line(file1,"\\\\")
latex_add_line(file1,"Correction for Measurement Error: & No & Yes \\\\")
latex_add_line(file1,"\\\\")
latex_add_line(file1,"\\textbf{Estimates of Voter Positions:}\\\\")
latex_table_entries(file1,rownames(tab2)[1:4],tab2[1:4,])
latex_add_line(file1,"\\\\")
latex_add_line(file1,"\\textbf{Estimates of Legislator Positions:}\\\\")
latex_table_entries(file1,rownames(tab2)[5:8],tab2[5:8,])
latex_add_line(file1,"\\\\")
latex_table_end(file1,caption="Validation of Estimates of the Positions of Voters and Legislators -- Estimates of the positions of the median voter in each state is correlated state ideology. Estimates of the positions of legislators are correlated with alternative estimates from \\citep{BattistaEtAl13} and \\citep{ShorMcCarty11}. We our estimates and the alternative estimates were commutated based on the same underlying data, we did not report correlations corrected for measurement error.",label="ValidationTable")

### TABLE 3 ###
bias.func <- function(x) { x[1] + x[2] * usmed - usmed }

# tax_govc and tax_govesta are only avaliable for 2 states, spend_agris in only avaliable for s
issues2 <- c("tax_alcohols","tax_capitalgainss","tax_cigarettes","tax_corps","tax_gass","tax_incomeLT75s","tax_incomeGT75s","tax_propertys","tax_saless","spend_enviros","spend_healthcares","spend_highereds","spend_k12s","spend_lawenforcements","spend_transinfras","spend_welfares")
tab3 <- matrix(rep("",2*(length(issues2)+3)*7),ncol=7)
for(i in 1:length(issues2))
{
	# correlation (corr. for m.e.)
	tab3[2*i-1,2] <- dig3(cor0(x,medvoter,subset=issue==issues2[i]))

	# correlation (corr. for m.e.)
	tab3[2*i-1,3] <- dig3(cor0(x,medvoter,x.se=x_se,y.se=medvoterse,subset=issue==issues2[i]))

	# responsiveness
	lm1 <- lm0(x ~ medvoter,subset=issue==issues2[i],w=1/x_se^2)
	if(lm1$N >= 30)
	{
		tab3[2*i-1,4] <- as.coef(lm1$coef[2],lm1$se[2])
		tab3[2*i,4] <- as.se(lm1$se[2])
	
		# bias	
		dm1 <- delta.method(bias.func,lm1$coef,lm1$V)
		tab3[2*i-1,5] <- as.coef(dm1$est,dm1$se)
		tab3[2*i,5] <- as.se(dm1$se)
	
		# RMSE
		tab3[2*i-1,7] <- dig3((mean0((x-medvoter)^2,subset=issue==issues2[i],w=1/x_se^2) - mean0(x_se^2,subset=issue==issues2[i],w=1/x_se^2) - mean0(medvoterse^2,subset=issue==issues2[i],w=1/x_se^2))^.5)
	} else {
		tab3[2*i-1,4] <- "n/a"
		tab3[2*i-1,5] <- "n/a"
		tab3[2*i-1,7] <- "n/a"
	}
	# N
	tab3[2*i-1,6] <- lm1$N
}
tab3[1,1] <- "Alcohol"
tab3[3,1] <- "Capital Gains"
tab3[5,1] <- "Cigarette"
tab3[7,1] <- "Corporate"
tab3[9,1] <- "Gasoline"
tab3[11,1] <- "Income, $<$ 75k"
tab3[13,1] <- "Income, $>$ 75k"
tab3[15,1] <- "Property"
tab3[17,1] <- "Sales"
tab3[19,1] <- "Environment"
tab3[21,1] <- "Health Care"
tab3[23,1] <- "Higher Education"
tab3[25,1] <- "K-12 Education"
tab3[27,1] <- "Law Enforcement"
tab3[29,1] <- "Transportation and Infrastructure"
tab3[31,1] <- "Welfare"
tab3[33,1] <- "\\textbf{Policy Index:}"
tab3[35,1] <- "\\textbf{Tax Index:}"
tab3[37,1] <- "\\textbf{Spending Index:}"

lm1 <- lm0(x_index ~ medvoter2,w=1/x_index_se^2)
tab3[33,4] <- as.coef(lm1$coef[2],lm1$se[2])
tab3[34,4] <- as.se(lm1$se[2])
tab3[33,2] <- dig3(cor0(x_index,medvoter2))
tab3[33,3] <- dig3(cor0(x_index,medvoter2,x.se=x_index_se,y.se=medvoterse2))
dm1 <- delta.method(bias.func,lm1$coef,lm1$V)
tab3[33,5] <- as.coef(dm1$est,dm1$se)
tab3[34,5] <- as.se(dm1$se)
tab3[33,6] <- lm1$N
tab3[33,7] <- dig3((mean0((x_index - medvoter2)^2) - mean0(x_index_se^2) - mean0(medvoterse2^2))^.5)

lm1 <- lm0(x_tax_index ~ medvoter2,w=1/x_tax_index_se^2)
tab3[35,4] <- as.coef(lm1$coef[2],lm1$se[2])
tab3[36,4] <- as.se(lm1$se[2])
tab3[35,2] <- dig3(cor0(x_tax_index,medvoter2))
tab3[35,3] <- dig3(cor0(x_tax_index,medvoter2,x.se=x_tax_index_se,y.se=medvoterse2))
dm1 <- delta.method(bias.func,lm1$coef,lm1$V)
tab3[35,5] <- as.coef(dm1$est,dm1$se)
tab3[36,5] <- as.se(dm1$se)
tab3[35,6] <- lm1$N
tab3[35,7] <- dig3((mean0((x_tax_index - medvoter2)^2) - mean0(x_tax_index_se^2) - mean0(medvoterse2^2))^.5)

lm1 <- lm0(x_spend_index ~ medvoter2,w=1/x_spend_index_se^2)
tab3[37,4] <- as.coef(lm1$coef[2],lm1$se[2])
tab3[38,4] <- as.se(lm1$se[2])
tab3[37,2] <- dig3(cor0(x_spend_index,medvoter2))
tab3[37,3] <- dig3(cor0(x_spend_index,medvoter2,x.se=x_spend_index_se,y.se=medvoterse))
dm1 <- delta.method(bias.func,lm1$coef,lm1$V)
tab3[37,5] <- as.coef(dm1$est,dm1$se)
tab3[38,5] <- as.se(dm1$se)
tab3[37,6] <- lm1$N
tab3[37,7] <- dig3((mean0((x_spend_index - medvoter2)^2) - mean0(x_spend_index_se^2) - mean0(medvoterse2^2))^.5)

latex_table_start(file1,5)
latex_add_line(file1," & \\textbf{Congruence} & \\textbf{Correlation} & \\textbf{Slope} & \\textbf{Bias} & \\textbf{N} \\\\")
latex_add_line(file1,"\\\\")
latex_add_line(file1,"\\textbf{Taxes:}\\\\")
latex_table_entries(file1,paste("\\hspace{.4cm}",tab3[1:18,1],sep=""),tab3[1:18,c(7,3:6)])
latex_add_line(file1,"\\textbf{Spending:}\\\\")
latex_table_entries(file1,paste("\\hspace{.4cm}",tab3[19:32,1],sep=""),tab3[19:32,c(7,3:6)])
latex_table_entries(file1,tab3[33:38,1],tab3[33:38,c(7,3:6)])
latex_table_end(file1,caption="Representation in Elections -- Congruence, Responsiveness, and Bias of Policy Outcomes to Voter Preferences. \\textit{Congruence} is the root mean square error between the estimated policy outcome or policy index and the median voter's position, corrected for measurement error. \\textit{Correlation} is the correlation between the estimated policy outcome or policy index and the median voter's position, corrected for measurement error. \\textit{Slope} is the slope of a linear regression where the dependent variable is a estimated policy outcome or policy index. \\textit{Bias} is the predicted value of from the regression evaluated at the position of the national median voter, minus the position of the national median voter. One star indicates statistical significance at the 5\\% level. Two stars indicates statistical significance at the 1\\% level. Three stars indicates statistical significance at the 0.1\\% level. A plus sign indicates statistical significance at the 10\\% level.",label="ResponsivenessTable")

### TABLE 4 ### 
lm1 <- lm0(medlegH2 ~ medvoter2,subset=state2!="AK"&state2!="HI"&state2!="NE"&state2!="US")
lm2 <- lm0(medlegS2 ~ medvoter2,subset=state2!="AK"&state2!="HI"&state2!="NE"&state2!="US")
lm3 <- lm0(gov2     ~ medvoter2,subset=state2!="AK"&state2!="HI"&state2!="NE"&state2!="US")

d1 <- delta.method(bias.func,lm1$coef,lm1$V)
d2 <- delta.method(bias.func,lm2$coef,lm2$V)
d3 <- delta.method(bias.func,lm3$coef,lm3$V)

tab4 <- models.to.latex(list("(1)"=lm1,"(2)"=lm2,"(3)"=lm3))
tab4 <- rbind(tab4[1:4,],rep("",3),tab4[5:6,],dig3(as.numeric(tab4[6,])^.5),rep("",3),dig3(c(d1$est,d2$est,d3$est)),c(as.se(d1$se),as.se(d2$se),as.se(d3$se)))
rownames(tab4) <- c("Constant","","Median Voter","","","N","$R^2$","Correlation","","Bias at U.S. Median","")

latex_table_start(file1,3)
latex_add_line(file1," & (1) & (2) & (3) \\\\")
latex_add_line(file1,"\\textbf{DV:} & \\textbf{House} & \\textbf{Senate} & \\textbf{Governor} \\\\")
latex_add_line(file1,"\\ & \\textbf{Median} & \\textbf{Median} \\\\")
latex_add_line(file1,"\\\\")
latex_table_entries(file1,rownames(tab4),tab4)
latex_table_end(file1,caption="Representation in Elections -- The dependent variables are the ideal points of the pivotal actors in each state and the independent variable is the ideal point of the state's median voter. One star indicates statistical significance at the 5\\% level. Two stars indicates statistical significance at the 1\\% level. Three stars indicates statistical significance at the 0.1\\% level. A plus sign indicates statistical significance at the 10\\% level.",label="ElectionRepresentationTable")


### TABLE 5 ###
lm1 <- lm0(x_index ~ st_ideo,subset=state2!="AK"&state2!="HI"&state2!="NE")
lm2 <- lm0(x_index ~ cont2,subset=state2!="AK"&state2!="HI"&state2!="NE")
lm3 <- lm0(x_index ~ contavg2,subset=state2!="AK"&state2!="HI"&state2!="NE")
lm4 <- lm0(x_index ~ contavg2,subset=state2!="AK"&state2!="HI"&state2!="NE"&!south2)
lm5 <- lm0(x_tax_index ~ contavg2,subset=state2!="AK"&state2!="HI"&state2!="NE"&!south2)
lm6 <- lm0(x_spend_index ~ contavg2,subset=state2!="AK"&state2!="HI"&state2!="NE"&!south2)
tab5 <- models.to.latex(list("(1)"=lm1,"(2)"=lm2,"(3)"=lm3,"(4)"=lm4,"(5)"=lm5,"(6)"=lm6))
tab5 <- rbind(tab5[1:8,],rep("",6),c("","","","X","X","X"),tab5[9:10,],dig3(as.numeric(tab5[10,])^.5))
rownames(tab5) <- c("Constant","","State Ideology Index","","Party Control (current)","","Party Control (20-year average)","","","\\textbf{South Excluded}","N","$R^2$","Correlation")

latex_table_start(file1,6)
latex_add_line(file1," & (1) & (2) & (3) & (4) & (5) & (6) \\\\")
latex_add_line(file1,"\\textbf{Policies:} & \\textbf{All}& \\textbf{All}& \\textbf{All}& \\textbf{All} & \\textbf{Tax} & \\textbf{Spending} \\\\")
latex_add_line(file1,"\\\\")
latex_table_entries(file1,rownames(tab5),tab5)
latex_table_end(file1,caption="State Government Composition and Policy Outcomes -- The dependent variables are policy indices. One star indicates statistical significance at the 5\\% level. Two stars indicates statistical significance at the 1\\% level. Three stars indicates statistical significance at the 0.1\\% level. A plus sign indicates statistical significance at the 10\\% level.",label="GovCompositionTable")

### TABLE 6 ### --- not based on calculations

### TABLE 7 ###

tab7 <- matrix(rep(NA,6*3),ncol=3)

tab7[1,2] <- dig2(mean0(ewm98,subset=state2!="US"))
tab7[2,2] <- dig2(mean0(ewm99,subset=state2!="US"))
tab7[3,2] <- dig2(mean0(ewm00,subset=state2!="US"))
tab7[4,2] <- dig2(mean0(ewm01,subset=state2!="US"))
tab7[5,2] <- dig2(mean0(ewm02,subset=state2!="US"))
tab7[6,2] <- dig2(mean0(ewm03,subset=state2!="US"))

tab7[1,3] <- dig2(sd0(ewm98,subset=state2!="US"))
tab7[2,3] <- dig2(sd0(ewm99,subset=state2!="US"))
tab7[3,3] <- dig2(sd0(ewm00,subset=state2!="US"))
tab7[4,3] <- dig2(sd0(ewm01,subset=state2!="US"))
tab7[5,3] <- dig2(sd0(ewm02,subset=state2!="US"))
tab7[6,3] <- dig2(sd0(ewm03,subset=state2!="US"))

tab7[1,1] <- as.per0(cor0(ewm00,ewm98,x.se=ewm00se,y.se=ewm98se,subset=state2!="US"),tex=TRUE)
tab7[2,1] <- as.per0(cor0(ewm00,ewm99,x.se=ewm00se,y.se=ewm99se,subset=state2!="US"),tex=TRUE)
tab7[4,1] <- as.per0(cor0(ewm00,ewm01,x.se=ewm00se,y.se=ewm01se,subset=state2!="US"),tex=TRUE)
tab7[5,1] <- as.per0(cor0(ewm00,ewm02,x.se=ewm00se,y.se=ewm02se,subset=state2!="US"),tex=TRUE)
tab7[6,1] <- as.per0(cor0(ewm00,ewm03,x.se=ewm00se,y.se=ewm03se,subset=state2!="US"),tex=TRUE)

latex_table_start(file1,3)
latex_add_line(file1,"\\textbf{Year} & \\textbf{Corr. w/ 2000} & \\textbf{Mean of} & \\textbf{Std. Dev. in} \\\\")
latex_add_line(file1,"& \\textbf{State Opinion} & \\textbf{State Opinion} & \\textbf{State Opinion} \\\\")
latex_table_entries(file1,1998:2003,tab7)
latex_table_end(file1,caption="Persistence of State Opinion",label="StateOpinionTable")

### TABLE 8 ###
tab8 <- cbind(HouseBridges,HouseRho,HousePerCorr,SenBridges,SenRho,SenPerCorr)
rownames(tab8) <- state2
latex_table_start(file1,6)
latex_add_line(file1," & \\multicolumn{3}{c}{\\textbf{House}} & \\multicolumn{3}{c}{\\textbf{Senate}} \\\\")
latex_add_line(file1," & Num. Bridges & $\\rho$ & Per. Corr.& Num. Bridges & $\\rho$ & Per. Corr. \\\\")
latex_add_line(file1,"\\\\")
latex_table_entries(file1,rownames(tab8),tab8)
latex_add_line(file1,"\\\\")
latex_table_end(file1,caption="Bridging the State Legislatures -- \\textit{Num. Bridges} refers to the number of legislators in the chamber that responded to the NPAT. $\\rho$ is the correlation between the NPAT ideal points and roll call ideal points for the chamber. \\textit{Per. Corr.} is the percent of roll call votes correctly predicted by a one-dimensional model in the chamber.",label="BridgeTable1")

### TABLE 9 ###
tab9 <- cbind(HouseNpatN,SenNpatN,GovNpatN)
rownames(tab9) <- state2
latex_table_start(file1,3)
latex_add_line(file1," & \\textbf{House} & \\textbf{Senate} & \\textbf{Governor} \\\\")
latex_add_line(file1,"\\\\")
latex_table_entries(file1,rownames(tab9),tab9)
latex_add_line(file1,"\\\\")
latex_table_end(file1,caption="NPAT Respondents --- The number of respondents from House candidates, Senate candidates, and gubernatorial candidates in each state. Responses come from the 1998 and 2000 NPATs in most states, with exceptions listed in Table \\ref{NPATYearsTable}",label="NPATTable1")

### EXTRA STATS ###
extra1 <- matrix(rep(NA,19*2),ncol=2)

# reliability of median voter estimate
extra1[1,1] <- "Reliability of Median Voter"
extra1[1,2] <- dig3(reli0(medvoter2,medvoterse2,subset=state2!="US"&state2!="AK"&state2!="HI"))

# libcon position on bridge votes
extra1[2,1] <- "Percent Liberal Responses, Democratic Candidates"
extra1[2,2] <- dig3(mean0(libcon,by=File3$Party,bys=c("D","d","r","R")))[1]

extra1[3,1] <- "Percent Liberal Responses, Democratic Voters"
extra1[3,2] <- dig3(mean0(libcon,by=File3$Party,bys=c("D","d","r","R")))[2]

extra1[4,1] <- "Percent Liberal Responses, Republican Voters"
extra1[4,2] <- dig3(mean0(libcon,by=File3$Party,bys=c("D","d","r","R")))[3]

extra1[5,1] <- "Percent Liberal Responses, Republican Candidates"
extra1[5,2] <- dig3(mean0(libcon,by=File3$Party,bys=c("D","d","r","R")))[4]

extra1[6,1] <- "Percent Liberal Responses, NAES Respondents"
extra1[6,2] <- dig3(mean0(libcon,subset=File3$Party=="d"|File3$Party=="i"|File3$Party=="r"))

# national median
extra1[7,1] <- "National Median Voter"
extra1[7,2] <- dig3(usmed)

# steny hoyer-joe biden distance
extra1[8,1] <- "Steny Hoyer-Joe Biden Distance"
extra1[8,2] <- dig3(abs(alpha4[fmatch("Cand_MMMUSS012_USS_1998",id4)] - alpha4[fmatch("Cand_MMMUSH096_USH_1998",id4)]))

# sussan collins-mitch mcconnell distange
extra1[9,1] <- "Susan Collis-Mitch McConnel Distance"
extra1[9,2] <- dig3(abs(alpha4[fmatch("Cand_MMMUSS029_USS_1998",id4)] - alpha4[fmatch("Cand_MMMUSS026_USS_1998",id4)]))

# steny hoyer-diane fienstein distange
extra1[10,1] <- "Steny Hoyer-Diane Feinstein Distance"
extra1[10,2] <- dig3(abs(alpha4[fmatch("Cand_MMMUSS007_USS_1998",id4)] - alpha4[fmatch("Cand_MMMUSH096_USH_1998",id4)]))

# EWM correlation between policy index and state opinion
extra1[11,1] <- "Erikson/Wright/McIver Policy-Opinion Correlation"
extra1[11,2] <- dig2(.68^.5)

# EWM rage of correlations
extra1[12,1] <- "Erikson/Wright/McIver Minimal Policy-Opinion Correlation"
extra1[12,2] <- dig2(min(c(.558^.5,.433^.5,.423^.5,.313^.5,.194^.5,.575^.5,.291^.5,.218^.5)))

extra1[13,1] <- "Erikson/Wright/McIver Maximal Policy-Opinion Correlation"
extra1[13,2] <- dig2(max(c(.558^.5,.433^.5,.423^.5,.313^.5,.194^.5,.575^.5,.291^.5,.218^.5)))

# Caughey-Warshaw comparisons
extra1[14,1] <- "Caughey-Warshaw Social Policy-Opinion Correlation"
extra1[14,2] <- dig2(cor0(cw_mass_soc,cw_pol_soc))

extra1[15,1] <- "Caughey-Warshaw Econ Policy-Opinion Correlation"
extra1[15,2] <- dig2(cor0(cw_mass_econ,cw_pol_econ))

# Erikson, wright, and McIver comparison
# on page 85, the report an R-squared of 0.218, suggesting the following correlation
extra1[16,1] <- "Erikson/Wright/McIver Tax Policy-Opinion Correlation"
extra1[16,2] <- dig2(.218^.5)

# grumbach comparison
extra1[17,1] <- "Grumbach Tax Policy-Opinion Correlation"
extra1[17,2] <- dig2(cor0(grum_tax,medvoter2))

# status quo vs. policy
extra1[18,1] <- "Status Quo-Policy Correlation"
extra1[18,2] <- dig2(cor0(x,s,subset=state!="US"))

# vote share in 2000 and 2004
extra1[19,1] <- "Right-wing Vote Share 2000-2004 Correlation"
extra1[19,2] <- dig2(cor0(rvs00,rvs04,subset=state2!="US"))

latex_table_start(file1,2)
latex_add_line(file1,"\\\\")
latex_table_entries(file1,extra1[,1],as.matrix(extra1[,2]))
latex_add_line(file1,"\\\\")
latex_table_end(file1)

latex_doc_end(file1)

